Technical  Report  SL-93-16 
September  1 993 


US  Army  Corps 
of  Engineers 

Waterways  Experiment 
Station 


AD-A271  540 

I  Mllti  ftl  MTU  iAMi  iiaa.  _ 

1 


Two-Dimensional  Finite  Element 
Analysis  of  Porous  Media 
at  Muitikiiobar  Stress  Levels 

by  Stephen  A.  Akers 
Structures  Laboratory 


%  'li 


Approved  For  Public  ReleeM:  Oietribution  It  UnHrrowd 


93  in 


057 


93-26102 

IIIIIIHII 


Prepared  for  Discretionary  Research  Program  • 

U.S.  Army  Engineer  Waterways  Experiment  Station 


The  contents  of  this  report  ire  not  to  be  used  for  advertising, 
publicition,  or  promotional  purposes.  Citation  of  trade  names 
does  not  constitute  an  ofTidal  endorsement  or  approval  of  the  use 
of  such  commercial  products. 


mwnDON  MCYOioMmi 


Technical  Report  SL-93-16 
September  1993 


TwO'Dimensional  Finite  Element 
Analysis  of  Porous  Media 
at  Multikllobar  Stress  Levels 

by  Stephen  A.  Akers 
Structures  Laboratory 

U.S.  Army  Corps  of  Engineers 
Waterways  Experiment  Station 
3909  Halls  Ferry  Road 
Vicksburg.  MS  39180-6199 


Final  report 

Aooreved  tor  oubiic  releossi  dvMbution  it  unlmilscl 


Prepared  for  Dbaetionary  Research  Program 

U.S.  Army  Engineer  Waterways  Experiment  Station 
3909  Halle  Ferry  Road.  Vicksburg.  MS  39180-6199 


WatioMyt  Exptrlimrrt  StMlon  Cataiogino^^ 

Akars.  Staphan  A. 

Two-dknanikmal  Mia  alamani  anatysis  o(  porous  madia  at  multidiobar  ttr^ 
lavais  /  by  Slaphan  A.  Akars :  praparad  tor  Otscralionary  Risaarch  Program. 

U.S.  Army  Engtoaar  Waiarways  Ei^wrimani  Station. 

105  p. :  M. ;  28  cm.  —  (Tachnicat  raport ;  SL-93-16) 

Indudas  btoiooraphical  ratoranoas. 

1.  Sols— Taalino.  2.  Porous  maiarials—Tastino.  S.Mattar— Pro(«artias— 
Mathamatical  modais.  4.  Finita  alamani  mathod.  I.  U.S.  Army  Enginaar  Mator- 
ways  Exparimani  Slattoix  H.  Diacrationary  Rasaaich  Program  (U.S.  Army 
Enginaar  Waiarways  Exparlmani  Station)  UlTWa.  IV.  Sariaa:  Tachnicd  raport 
(U.S.  Army  Enginaar  Watamvays  Exparimant  Station) ;  SL-93-16. 
TA7W34no.SL-03>16 


Contents 


Preface  .  v 

1 - Introduction  .  I 

Background .  1 

Approach .  2 

Purpose  and  Scope  .  3 

2- Finite  Element  Model  .  4 

Introduction  .  4 

Background .  4 

Finite  Element  Fommlation .  6 

Effective  stresses  and  strains .  6 

Finite  element  equations . 7 

Equations  for  residual  forces . II 

Constitutive  Models . 15 

Element  Implemented  into  JAM . IS 

Sununary  . 16 

3- The  Cap  Model  and  Its  ImplementaticMi . 17 

Introduction  . 17 

Background . 17 

Loading  Functions  and  Row  Rule . 19 

Derivation  of  Incremental 

Elastic-Plastic  Stress-Strain  Relations  . 21 

Elastic-Plastic  Constitutive  Matrix . 25 

Plastic  Hardening  Modulus . 27 

Cap  Model  In^lemented  into  JAM  . 28 

Numerical  Implementation  of  the  Cap  Model . 33 

Itttroduction  . 33 

Elastic  algorithm  . 34 

Failure  envelope  algorithm  . 34 

Cq}  algCNrithm . 40 

Tensile  algorithm . 42 

Implementatkm  of  Cap  Model  into  JAM . 43 

Verification . 44 


Summary 


49 


4- Equations  of  State  for  Air, 

Water,  and  Solids  . 50 

Introduction  . 50 

Equation  of  State  for  Water . 50 

Air-Water  Compressibility  . 52 

Background . 52 

Boyle’s  and  Henry’s  laws . 53 

Derivation  of  equations  . 54 

Equation  of  State  for  Solids . 59 

Summary  . 62 

5- Features  and  Verification  of  FE  Program . 63 

Introduction  . 63 

Additional  Features  of  FE  Program . 63 

Restart  feature . 63 

Postprocessing . 64 

Plane  and  Axisymmetric  Verification 

Problems . 64 

Consolidation  Problems  . 70 

Cryer  Problem .  73 

Summary  .  75 

6- Numerical  Simulations .  76 

Introduction  .  76 

Salem  Limestone  .  76 

Simulations .  76 

Process  .  76 

Drained  limestone  simulations  .  77 

Undrained  limestone  simulations .  79 

Test  Specimen  Sinailations  .  84 

FE  grid .  84 

SinuUation  of  drained  triaxial 

compression  test  .  84 

Simulation  of  consolidated  undrained 

triaxial  compression  test  .  88 

Summary  . 91 

7- Summary .  92 


References 


93 


IV 


Preface 


This  investigation  was  sponsored  by  the  U.S.  Army  Engineer  Waterways 
Experiment  Station  (WES)  under  the  Discretionary  Research  Program. 

The  work  was  conducted  by  Mr.  Stefriien  A.  Akers,  Geomechanics 
Division  (GD),  Structures  Laboratory  (SL).  WES,  during  the  period 
January  1991  through  March  1993,  under  the  general  direction  of  Dr.  J.  G. 
Jackson,  Jr.,  Chief,  GD,  and  Mr.  J.  Q.  Ehrgott,  Chief,  Operations  Group, 
GD.  Mr.  Bryant  Mather  is  Director,  SL. 

At  the  time  of  publication  of  this  report.  Dr.  Robert  W.  Whalin  was  the 
Director  of  WES.  Commander  was  COL  Bruce  K.  Howard,  EN. 


1 


Introduction 


Background 

The  personnel  of  the  Geomechanics  Division,  Structures  Laboratory.  U  S. 
Army  ^ineer  Waterways  Experiment  Station  (WES)  are  responsible  for 
research  and  development  in  the  general  field  of  soil  and  rock  dynamics.  Our 
primary  interest  is  in  the  response  of  earth  and  eaith-stnumire  systems 
subjttted  to  intense  transient  (blast-type)  loadings.  An  analysis  of  these 
systems  is  typically  conducted  in  three  different  phases.  First,  laboratory  tests 
:irc  conduct^  on  the  geologic  materials  of  interest  in  order  to  develop  a  data 
base  of  c(Mq>osition  and  mechanical  properties;  then,  based  upon  this  data 
base,  a  set  of  recommended  material  properties  is  developed  for  the  consti- 
nttive  modelen.  In  the  second  frfuse.  the  modelers  fit  one  or  more  consti¬ 
tutive  models  to  the  recommended  material  properties.  Last,  fmite  element 
(FE)  or  fmite  difference  codes  are  used  by  the  ground  shock  calculators  to 
analyze  the  responses  of  these  systems. 

The  Geomechanics  Division  is  frequently  asked  to  conduct  mechanical 
property  investigations.  In  performing  these  investigations,  we  have  tested 
and  characterized  many  types  of  materials.  These  materials  generally  fall  into 
the  following  groups;  naoisi  and  fiilly  saturated  cohesionleu  soils,  daen 
alluvituns,  natural  and  remolded  clays,  clay  diales.  soil  and  rock  "matching'' 
grouu,  and  a  variety  of  competent  rocks.  As  basing  and  attack  scenarios  of 
the  Department  of  Defense  become  more  elaborate,  and  as  the  analysis  tech¬ 
niques  of  modelen  and  ground  shock  calculaton  become  more  refined, 
greater  demands  are  placed  on  the  engineer  who  is  asked  to  perform  and 
analyze  the  laboratory  mechanical  property  tesu  in  order  to  provide  recom¬ 
mended  material  properties.  Modelen  and  calculaton  are  now  requesting 
tmal-stress  mKhanical  dau  at  stress  levels  of  several  kiltrfian. 

Coo^licated  stress-  and  strainipath  tests  are  frequently  included  in  their  lisu 
of  dnired  muerial  response  tesu.  Greater  erof^is  in  effeaive-stress 
material  properties  is  now  evident  in  the  ground  shock  community. 

Due  to  tlK  unconventional  nature  of  many  of  the  requested  tesu.  the 
engineer  performir^  and  analyzing  the  tesu  is  sometimes  uncertain  aboM 
existing  l^ioratory  equipment,  i.e..  whether  it  restricu  the  types  of  tesu  that 
can  be  conducted.  The  engineer  may  question  the  measured  laboratory 
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responses;  are  they  theoretically  realistic?  An  engineer  may  have  specific 
questions  such  as:  (a)  what  effect  will  small  amounts  of  air-filled  porosity 
have  on  material  properties,  (b)  what  loading  rates  are  appropriate  for  con¬ 
ducting  truly  drained  tests  or  undrained  tests  with  meaningful  pore  pressure 
measurements,  and  (c)  how  does  one  calculate  effective  stresses  at  kilobar 
stress  levels  for  rock-like  materials?  In  some  situations,  an  engineer  respon¬ 
sible  for  recommending  material  properties  may  only  have  low  pressure  (less 
than  a  kilobar)  total-stress  and  effective-stress  data  ftom  which  to  extrapolate 
multikilobar  material  responses. 

An  engineer  would  have  a  tremendous  advantage  if  a  numerical  tool  were 
available  with  which  to  verify  laboratory  test  results  or  to  predict  unavailable 
laboratory  test  data.  The  appropriate  numerical  tool  should  give  an  engineer 
the  capability  to  calculate  both  total-  and  effective-stress  material  responses. 
This  numerical  tool  would; 

1 .  calculate  strains,  total  and  effective  stresses,  and  pore  fluid  pressures 
for  fiilly-  and  partially-saturated  porous  media, 

2.  calculate  the  time  dependent  flow  of  pore  fluids  in  porous  media, 

3.  model  nonlinear  irreversible  stress-strain  behavior,  including  coupled 
shear-induced  volume  change,  and 

4.  simulate  the  effect  of  nonlinear  pore  fluid  compressibility  and  the 
contribution  of  the  compressibility  of  the  grain  solids  for  stresses  up  to 
several  kilobars. 

The  FE  code  JAM  incorporates  all  of  the  above  features.  The  code  simu¬ 
lates  quasi-sutic,  axisymmetric,  laboratory  mechanical  property  tests,  i.e.,  the 
laboratory  tests  are  analyzed  as  boundary  value  problems.  Features  I  and  2 
were  incorporated  into  the  code  using  modified  formulations  of  Biot's  coupled 
theory  as  advanced  by  investigators  such  as  Zienkiewicz  (1985a)  and  Lewis 
and  Schrefler  (1987).  An  elastic-plastic  strain-hardening  cap  model  calculates 
the  time-independent  skeletal  responses  of  the  porous  solids.  This  enables  the 
code  to  model  nonlirtear  irreversible  stress-strain  behavior  and  shear-induced 
volume  changes.  In  order  to  accurately  model  the  total-  arxi  effective-stress 
responses  of  multikilobar  laboratory  tests,  fluid  and  solid  compics sibilities 
were  incorporated  into  the  code.  Following  the  concept  used  by  Chang  and 
Duncan  (1983),  partially-saturated  materials  were  simulated  with "hontoge- 
nized*  compressible  pore  fluid. 


Approach 

To  develop  the  FE  code,  four  tiujor  tasks  were  conqileted.  They  were. 

I )  a  cap  model  was  incorporated  iiHo  an  exuiing  (wo-dimensional  finite 
elemeiM  code  PLAST  and  numerical  experiments  were  conduaed  to 
verify  iu  implemeittation; 
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2)  a  modified  version  of  Biot’s  coupled  theory  was  implemented  into  the 
code  produced  in  step  1  above; 

3)  a  data  base  of  drained  and  undrained  mechanical  properties  was  ob¬ 
tained  for  a  limestone  with  a  porosity  of  13.S%;  it  included  me:hani- 
cal  properties  from  tests  conducted  at  stress  levels  of  up  to  6  kilobars 
and  recommended  properties; 

4)  using  the  recommended  properties  from  step  3,  a  cap  model  was  fit  to 
the  drained  properties;  the  undrained  stress-strain  and  pore  pressure 
responses  of  the  material  were  calculated  with  the  FE  code  and  then 
compared  to  the  undrained  material  responses. 


Purpose  and  Scope 


The  purpose  of  this  report  is  to  document  the  features  and  algorithms 
implemented  into  the  FE  code  JAM.  Chapter  2  describes  the  FE  model  im¬ 
plemented  into  JAM  and  briefly  documents  the  constitutive  models  available  in 
the  code.  The  essential  features  of  the  cap  model  are  reviewed  and  the  steps 
required  to  implement  the  cap  model  into  the  FE  code  JAM  are  summarize 
in  Chapter  3.  The  equations  of  state  for  air.  water,  and  grain  solids  are  docu¬ 
mented  in  Chapter  4.  and  the  equations  for  compressibility  of  an  air-water 
mixture  are  developed.  Chapter  S  describes  features  in  the  FE  program  not 
introduced  in  earlier  chapters  and  presents  solutions  from  sever^  verification 
problems  as  proof  that  the  program  works  correctly.  Numerical  simulations 
of  limestone  behavior  under  drained  and  undrained  boundary  conditions  are 
presented  in  Chapter  6.  The  final  chiqKer  summarizes  the  results  of  this 
research  effort. 
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2  Finite  Eiement  Modei 


Introduction 

This  ctu|)ter  describes  the  FE  model  implemented  into  JAM  and  briefly 
documents  the  constitutive  models  available  in  the  code.  The  work  of  Biot 
and  other  investigators  is  described,  followed  by  a  discussion  of  the  equations 
implemented  by  Lewis  and  Schrefler  and  modifications  that  must  be  made  to 
those  equations.  In  addition,  the  equations  are  derived  for  the  residuals. 
Finally,  the  five  constitutive  models  available  in  JAM  for  modelling  skeletal 
behavior  are  described. 


Background 

In  1941.  Biot  published  his  three-dimensional  theory  of  consolidation  for 
static  loading.  In  his  theory.  Biot  coupled  the  solution  of  the  equations  of 
pore  fluid  diffusion  with  the  equations  of  defonnati(m  for  the  porous  solids. 

He  was  thus  able  to  calculate  time-dependent  displacemenu,  strains,  pore  fluid 
pressures,  and  effeaive  stresses.  Biot  made  the  following  assumptions  in  his 
fonnulation:  (1)  the  material  was  isotropic.  (2)  the  material  was  linear  elastic. 
(3)  small  strains  were  applicable.  (4)  the  pore  water  was  incompressible.  (S) 
the  pore  water  could  contain  air  bubbles,  and  (6)  flow  of  the  pore  water  (4)ey- 
ed  Darcy's  Law.  In  subsequent  papers.  Biot  extended  his  theory  to  inclutte 
anisotropic  materials,  viscoelastic  materials,  and  dynamic  processes  (Biot 
1955;  1962). 

With  the  rapid  development  of  digital  con^uters  and  advances  in  numerical 
techniques  such  as  the  finite  element  method,  nuuiy  investigators  etqianded 
Biot’s  theory  in  attempts  to  model  more  realistic  and  nx>re  complex  problems. 
Sandhu  and  Wilson  (1969)  were  the  firtt  to  use  finite  donem  techniques  with 
Bkx’s  original  formulation  to  solve  initial  boundary  value  problents.  They 
applied  variatitmal  principles  to  the  field  equations  of  fluid  flow  in  a  fully 
saturated  porous  elastic  continuum,  and  that  used  the  finite  element  met^  to 
numerically  solve  the  resulting  coupled  equati<ms. 
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Ghaboussi  and  Wilson  (1972,  1973)  developed  a  variational  formulation  of 
Biot’s  dynamic  field  equations  for  saturated  porous  elastic  solids.  Their  finite 
element  formulation  allowed  for  the  compressibilities  of  both  the  fluid  and 
solid  phases.  Their  methods  were  applicable  to  dynamic  soil -structure  inter¬ 
action  and  wave  propagation  problems  in  saturated  geologic  media. 

Zienkiewicz  and  his  colleagues  have  written  extensively  about  their  use  of 
modified  versions  of  Biot’s  formulation  to  solve  consolidation,  liquefaction, 
and  wave  propagation  problems  in  fluid  saturated  porous  materials  (Simon  et 
al.  1986a  and  1986b;  Zienkiewicz  198Sa;  Zienkiewicz  et  al.  1980; 
Zienkiewicz  and  Shiomi  1984).  They  have  incorporated  several  different  non¬ 
linear  constitutive  models  into  their  numerical  codes.  For  example,  Lewis  et 
al.  (1976),  used  a  hyperbolic  constimtive  relationship  to  model  the  skeletal 
response  of  the  solids.  They  incorporated  fluid  and  solid  compressibilities, 
creep,  and  void  ’atio  dependent  permeability  into  their  code.  Zienkiewicz  and 
his  colleagues  have  also  demonstrated  a  c^ability  to  solve  dynamic  problems 
such  as  ground  motions  due  to  earthquakes.  For  example,  Zienkiewicz  et  al. 
(1982),  modified  the  Critical  State  m^el  to  include  a  Coulomb-type  failure 
surface  and  incorporated  a  "cumulative  densification*  feature  into  the  consti¬ 
tutive  model  that  allowed  pore  fluid  pressures  to  increase  with  increasing 
numbers  of  load-unload  cycles.  They  then  demonstrated  the  utility  of  their 
approach  when  they  used  their  code  to  approximate  the  earthquake  induced 
displacements  and  pore  pressures  within  the  Lower  San  Fenu^o  Dam. 

Other  investigators  have  also  expanded  Biot’s  formulation  with  nonlinear 
constitutive  models.  Oka  et  al.  (1986)  developed  an  elasto-viscoplastic  consti¬ 
tutive  model  for  clay  and  used  Biot’s  theory  to  study  the  two-dimensional  con¬ 
solidation  response  of  sensitive  and  aged  clay  deposits.  They  demonstrated 
the  capability  of  their  code  to  simulate  the  consolidation  response  of  clay 
deposits  during  die  construction  phase  of  embankments.  Chang  and  Duncan 
(1983)  used  a  modified  Cam  Clay  constitutive  model  in  their  analyses  of  earth 
Structures  constructed  of  con^cted,  partially  sanitated  clay  soils.  They  also 
applied  the  concept  of  a  "homogenized  pore  fluid*  to  Biot’s  formulation  in 
order  to  model  partially  saturated  clay  soils.  With  this  implementation,  they 
were  able  to  treat  partially  sanitated  soils  as  two-phase  materials,  i.e.,  solids 
and  compressible  pore  fluid,  instead  of  using  a  more  theoretically  rigorous 
approach  involving  three-phase  materials. 

Lewis  and  Schrefler  (1987)  extetxled  Biot’s  formulation  to  include  the 
governing  (»]uations  for  single  phase,  multiphase,  and  saturated-unsamrated 
flow  in  a  deforming  porous  solid.  They  discussed  finite  element  procedures 
for  both  the  space  and  itme  discretization  aspects  of  consolidation  problems. 
They  also  present«i  linear  elastic  and  nonlinear  constitutive  relationships;  the 
nonlinear  models  included  the  hyperbolic  model  and  incremenul  elastic-plastic 
models  such  as  the  Critical  State  models. 
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Finite  Element  Formulation 


Effvctivt  str«ts«s  and  strains 

For  a  nonlinear  material  not  susceptible  to  creep  strains,  a  general  stress- 
strain  relation  can  be  written  as 

2. 

“  ^ijkl 


where  a^j'  is  the  matrix  of  effective  stresses,  is  the  tangential  stiflness 

matrix  or  constitutive  matrix,  is  the  matrix  of  total  strains,  is  a 
matrix  of  strains  due  to  the  compression  of  the  grains  by  the  pore  fluid  and 

dt is  a  matrix  of  effective  strains.  The  matrix  dt h  is  evaluated  as: 


2.2 


where  is  the  bulk  modulus  of  the  grains,  r  is  the  pore  fluid  pressure,  6^  is 
the  Kronecker  delta  defined  by 


1  if  i j 
0  if  i  j 


and  an  engineering  mechanics  sign  convontion  is  used  in  which  compression  is 
negative. 

The  purpose  of  the  term  dth  is  illustrated  in  the  following  example.  If  a 
porous  specimen  surrounded  by  a  pervious  membrane  was  placed  into  a  pres¬ 
sure  vessel  and  a  pressure  of  several  hundred  megapascals  was  applied,  a 
volume  decrease  would  be  measured  in  the  specimen  due  to  the  compression 
of  the  grains  However,  since  the  total  strain  dt  is  equal  to  the  volume 

strain  due  to  grain  compressibility,  i.e..  dt °  rhe  effective  strains  and 
therefore  the  effective  stresses  within  the  specimen  are  zero.  With  the  term 

dt  h  included  in  the  general  stress-strain  relation,  effective  stiesses  can 
singly  be  calculated  as; 

‘  2.3 

Under  drained  boundary  conditions,  the  toul  and  effective  strains  are  «]ual. 
and  the  total  and  effective  stresses  are  equal 
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In  Equation  2.3,  no  factor  need  be  applied  to  the  pore  pressure  term  to 
account  for  grain  compressibility,  which  was  a  method  proposed  by  Skempton 
(1960).  After  2q>prq)riate  manipulation,  the  above  equations  will  yield 
Skempton' s  equation 


Ap'  =  Ap  - 


2.4 


where  p  is  pressure,  AT  and  are  the  bulk  modultis  of  the  skeleton  ami  grain 
solids,  respectively. 


Rnlt*  tlcment  •quations 

The  general  equations  developed  from  the  spacial  discretization  of  the 
equilibrium  and  continuity  equations  have  been  documented  by  a  large  number 
of  investigators  (Lewis  and  Schrefler  1987;  Zienkiewicz  198Sa).  Lewis  and 
Schrefler  devel(^>ed  the  following  equations: 

(/:]«  -  (f-l  »  -  -  0  2.5 


and 

\H]r  -  (5]ir  -  [ifii  -  Q-0 
where  is  the  tangent  stiffness  matrix  of  the  solid  phase. 


2.6 


K  -  (  B^Dj-BdQ 
0 


2.7 


(L]  is  the  coupling  matrix  between  the  solid  and  fluid  phases. 


L  »  [  B^mNdQ  -  [ 
•'a  •'a 


3K. 


NdQ 


2.8 


[H]  is  the  permeability  matrix  of  the  porous  skeleton. 


«.  j 


(VN)^  i  VNdQ 


2.9 


\S\  is  the  con^ressibility  matrix. 


S  »  [  N^s  N  dQ 
Q 


2.10 
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in  which  the  scalar  s  is  evaluated  as 


=  ^  +  A 


1 


Drin 


2.11 


«/  i3Kgr 

R  is  the  external  force  vector. 


/?  =  [ 


^dbdQ*  [  N^dtdT 
Q  '^r 


2.12 


Q  is  the  vector  of  boundary  flows. 


V pgh  dfi 


2.13 


and  the  superinqxjsed  dot  indicates  a  time  derivative.  In  the  above  equations, 

{(  is  the  vector  of  displacement  increments,  it  is  the  vector  of  pore  pressure 
increments,  D  is  the  elastic-plastic  constitutive  matrix,  m  is  the  matrix 
equivalent  of  the  Kronecker  delta,  is  the  strain-displacement  matrix,  N  is  the 
matri.x  of  displacement  shape  functions,  N  is  the  matrix  of  pore  pressure 

shape  functions,  6  is  a  vector  of  body  forces,  /  is  a  vector  of  surface  trac¬ 
tions,  k  is  the  absolute  permeability  matrix  of  the  material,  (i  is  ihe  dynamic 
viscosity  of  the  pore  fluid,  and  are  the  bulk  modulus  of  the  grain 
solids  and  pore  fluid,  respeaively,  ^  is  the  porosity  of  the  material,  q  is  the 
vector  of  applied  fluid  flux,  and  p ,  g  and  h  are  fluid  density,  gravity,  and 
elevation  head,  respectively. 

Lewis  and  Schrefler  (1987)  describe  in  detail  the  components  that  comrib- 
ute  to  the  rate  of  fluid  accumulation.  However,  one  of  the  terms  in  their 
formulation  is  misleading  if  not  incorrect.  The  following  discussion  expands 
on  their  analysis  and  then  shows  why  their  formulation  requires  modification. 
For  the  material  and  conditions  of  interest,  Lewis  and  Schrefler  specify  that 
four  volumetric  strain  con^nents  must  be  evaluated;  they  are  the  volume 

strain  of  the  porous  matrix  volume  strains  of  the  grain  solids  < f  and 

the  pore  fluid  c  /,  and  a  component  of  solid  volume  strain  due  to  the 
applied  effective  stresses.  To  simplify  the  discussion  of  fluid  accumulation, 
consider  i)m  following  test  conditions.  A  fully  sanirated  porous  material  is 
contained  in  a  sample  chamber,  which  has  frictionless  sides,  and  is  loaded  by 
a  frictionless  piston  with  area  A.  A  flow  meter  attached  to  the  piston  indicates 
the  direction  of  fluid  flow  into  (positive  flow)  or  out  of  (negative  flow)  the 
sample.  The  tq)  surface  of  the  piston  is  loaded  by  a  chamber  pr»sure  P|  and 
a  second  fluid  pressure  P2  is  applied  through  the  flow  meter  to  the  pore  fluid 
within  the  san^le.  The  fluid  pressures  P|  and  P2  control  the  total  and  effec¬ 
tive  stresses  within  the  san^le.  Recall  that  for  a  porous  material  with  a 
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volume  of  unity,  the  volume  of  the  voids  is  n,  and  the  volume  of  the  solids  is 
1-n. 


If  P,  is  inaeased  and  P2  held  constant,  the  effective  stress  in  the  specimen 
will  increase,  the  material  will  compress,  and  pore  fluid  will  flow  out  of  the 
specimen.  Since  the  specimen  is  fully  samrated  with  fluid,  the  rate  of  change 
in  fluid  accumulation  is  equal  to  the  volumetric  strain  of  the  porous  matrix 
and  may  be  written  as; 


^  g  2.14 

dt  ~  dt  ~  dt 


If  P]  and  P2  are  increased  at  the  same  rate,  the  effective  stress  within  the 
specimen  will  not  change.  However,  both  the  grain  solids  and  the  pore  fluid 
will  compress  due  to  the  change  in  pore  fluid  pressure.  For  these  conditions, 
fluid  will  flow  into  the  specimen.  Let  us  evaluate  the  solid  and  fluid  compo¬ 
nents  separately.  The  volui  strain  in  the  solids  due  to  an  applied  pressure  is 
expressed  as: 


2.15 


where  is  the  bulk  modulus  of  the  grains  and  is  the  volume  of  the  sol¬ 
ids.  From  this  equation,  we  can  express  the  volume  change  within  the  unit 
volume  due  to  the  compression  of  the  grains  as; 

d(l  ^  -n)d(,  ^  (l-n)  2.16 

where  dP  has  been  replaced  by  the  pore  fluid  pressure  dir. 

In  the  same  manner,  the  volume  change  within  the  unit  volume  due  to  the 
compression  of  the  pore  fluid  is  written  as. 


where  Kj  is  the  bulk  modulus  of  the  pore  fluid. 

Lewis  and  Schrefler  also  evaluate  the  component  of  volume  strain  due  to 
the  compression  of  the  grains  caused  by  the  increase  in  effective  stress.  A 
similar  analysis  was  developed  by  Bishop  ( 1973).  Refer  again  to  the  test 
configuration  and  the  example  in  which  P|  was  increased  and  P2  held 
constant  For  a  sutistically  random  distribution  of  pore  space  within  the 
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specimen,  the  area  of  solids  (A,)  on  any  plane  through  the  specimen  will  be 

A^^{\-n)A  2.18 


where  n  is  the  porosity  of  the  specimen  as  previously  defined  If  do is  the 
average  normal  effective  stress  on  any  surface,  then  da I  {  \  -  n)  is  the 
average  normal  stress  in  the  solids.  Using  the  concept  implied  by  Equa¬ 
tion  2.  IS,  the  pressure  applied  to  the  solids  can  be  expressed  as: 

dP  -  Si,  .  —IL  219 

3(l-n)  3(l-n) 


which  when  substituted  into  Equation  2.16  gives; 


da./ 


^  V,  (l-n)V  3A:^(l-n) 

The  component  of  volume  strain  due  to  the  applied  elective  stresses  for  the 
unit  volume  is  then: 

...  _  dv  .  <'<’«■  ,  , 


The  components  of  the  volume  strain  are  now: 


and  the  expressions  for  each  component  when  substituted  gives; 

3  (oul  «  matrix  a  •  f 

°***  ^  _  1  ^  I  -  n  ^  n  dr  2.23 

"IT  ‘  dr  “  3Kg  dt  ~K^  Tf  If 


The  term  in  dispute  is  the  conqx)nent  of  volume  strain  due  to  the  applied 
effective  stresses.  Grain  compression  due  to  changes  in  effective  pressure  is 
already  incorporated  into  all  skeletal  constiniiive  models  by  default.  In  an 
analysis  of  drained  test  dau.  the  skeletal  strains  are  not  decoupled  from  the 
grain  strains,  and  the  sum  of  the  two  is  always  measured  by  strain  gauges  or 
deformeters  Therefore,  both  components  are  included  when  a  constitutive 
model  is  fit  to  drained  hydrostatic  loading  data.  A  similar  argument  can  be 
made  by  examining  Equation  2.23.  Duhng  a  drained  test.  Equation  2.2? 
would  reduce  to  the  following 
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Clearly,  this  is  in  error;  the  only  strains  iiKluded  in  the  above  expression 
should  be  the  skeletal  strains.  Therefore,  one  should  only  include  grain  com¬ 
pression  due  to  changing  pore  pressures  in  the  final  formulation.  The  neces¬ 
sary  modifications  were  made  to  Equations  2.8  and  2.11,  which  arc  rewritten 
below 


and 


L  =  [  B'^mNdQ 
•*0 


2.25 


s 


1 


K 


g 


2.26 


Equations  for  residual  forces 

Although  numerous  papers  pertaining  to  the  FE  equations  governing  pore 
fluid  flow  in  a  deforming  porous  solid  are  available,  none  outline  the  equa¬ 
tions  required  to  calculate  the  residual  forces.  In  this  section,  these  equations 
are  developed  for  a  nonlinear  incremental  finite  element  program  that  employs 
a  modified  Newton-Raphson  solution  scheme. 

The  time  integration  of  Equations  2.5  and  2.6  is  performed  using  the  fol¬ 
lowing  approximation; 

X  dt  a  A/  '*^x  *  ( 1  "  «)  A/  \  2.27 


for  0  ^  a  ^  1.  From  Equation  2.27,  the  following  are  developed; 


t*a£u . 


A/ 


2.28 


and 


I  •aAf 


X  =  ( I  -  o)  'x 


a 


f ‘Af, 


2.29 
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Table  2.1  gives  the  most  common  difference  schemes  adopted  by  the  selec¬ 
tion  of  a  given  value  of  a.  Equation  2.5  may  be  written  for  a  given  time 

Table  2.1. 

Time  Integration  Parameters 


ValiM  a 


Diffaranc*  Schama 


StabOitv 


Forward  or  Euler 


Crank-Nicolson 


G  alar  kin 


Backward 


Conditionally 

Unconditionally 

Unconditionally 


Conditionally 


t+aAt  as: 


r  +  aA/^  _  l*aAl^  _  t*aiilp 


[L]  . 


Equations  2.28  and  2.29  are  introduced  into  Equation  2.30  to  produce  the 
following 


-  ‘u 


I  I 

T  -  't 


'*^R-  'R 


Multiplying  by  A/  and  collecting  terms,  one  obuins 


,«A/[^]| /♦A/^  _  ,*A/^  _ 


,*AiR_  Iff 


In  general,  Equation  2.32  represents  nonlinear  behavior.  The  relationship 
may  be  linearized  with  the  following  expressions: 

/«A/y(()  .  /•A/yii-i)  ^  gyd)  2 


/♦A/^(i)  ^  /♦A/y(i-l)  ^  5y(') 


where  i  represents  the  current  iteration,  and  the  initial  conditions  are 

r*A/j^(0)  _  iji  jmd  /•A/jj.(0)  .  1^  This  linearization  can  be  used  as  the 
first  step  in  a  Newton-Raphson  iteration  (Bathe  1982)  If  Equations  2.33 


12 


Chcpicf  2  FirMte  Element  Model 


and  2.34  are  substituted  into  Equation  2.32  and  the  terms  for  iterations  (i-l) 
and  (0)  are  brought  to  the  right  hand  side  of  the  equation,  then  one  obtains: 

-  ‘R 

-  ] I -  '♦A/j^(0)|  2.3 

^  f+Ar  j  I /+A/^(i-l)  _  »+A/^(0)| 


Recognizing  that 


Equation  2.3S  may  be  written  in  terms  of  the  incremental  or  accumulative 
stresses  and  pore  pressures  as 


1)  ^  l♦A/^(l-l) 

where  ''^dV  and 

Equation  2.37  is  one  of  the  two  equations  required  to  solve  for  pore  fluid  flow 
in  a  deforming  nonlinear  porous  solid.  The  second  equation  is  developed 
from  Equation  2.6  and  is  written  at  time  t+aAf  as: 


-  [LI 


jT  l*a^^  ^  t*aAtQ 


•  • 


Equations  2.28  and  2.29  are  introduced  into  Equation  2.38  to  produce 


[f/]{(l-a)  'x  >  a  "^x}  -  [5] 

at 


ifl 


(I  -a)  'Q  *  a 
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Collecting  terms  and  multiplying  by  A/  one  obtains; 

{-[S]  +  aAr[Hl}  ‘*^‘u-'u) 


+  {[5]  +  (l-a)Ar[Wl}  'ir  =  P 


where  P  =  Ar{(l  -a)  ‘Q  +  a  If  the  following  equations 

5  =  -[5]  ♦  aAt[H] 


H  ^  [S]  *  (l-a)At[ff]  =  At[H]  -  S 

are  substituted  into  Equation  2.40  for  simplification,  then  one  obtains 
S  -  [Lf  '*^u  =  P  -  H  ‘r  -  [Lf  ‘u  . 


Equation  2.43  must  now  be  formulated  for  general  nonlinear  behavior 
using  the  same  linearization  process  applied  to  Equation  2.32,  i.e..  Equa¬ 
tions  2.33  and  2.34  must  be  substitute.  This  operation  produces 


_  I*  ♦Aiy(O)  _  i+A/i^jr  i+A/j^co) 


Collecting  the  6u  and  6t  terms  on  the  left  hand  side  of  the  equation,  one  gets 

2.45 

(i-l)  ^ 

Equations  2.41  and  2.42  can  be  used  to  eliminate  the  terms  ^  and  5  from  the 
right  hand  side  of  Equation  2.45  to  produce 

2.46 

Where  '-^G*'  **  «  {  -  ♦  aA/ //f ' ‘>1 

jr(i  -  1)  ^♦A/^J^(l-I) 

/♦A/^^(«-l)  .  /♦A/^(»-l)  „  I'Al^tO) 

fA/Ay('-l>  .  fAlj^(«-l)  .  fAlj^(O) 

».aia„(0)  .  0  and  »  0  . 
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Equations  2.37  and  2.46  are  written  in  matrix  form  for  increment  r+Ar  as: 


K  -  L 

-  -  S  *  aAtH 


5  »('■) 


jf  _ 


2.47 

This  is  the  system  of  equations  that  must  be  solved  to  calculate  displacements 
and  pore  fluid  pressures  in  a  deforming  porous  solid. 


Constitutive  Models 

Four  of  the  five  constitutive  models  available  in  the  FE  code  JAM  were 
deve  oped  by  Owen  and  Hinton  (1980).  These  models  were  implemented  by 
Owen  and  Hinton  in  the  FE  code  PLAST,  which  was  the  original  FE  code  on 
which  JAM  was  built.  Each  of  the  four  models  were  inq)lemented  as  elastic- 
plastic  models  with  linear  strain  hardening,  and  each  has  a  different  yield 
criteria.  The  four  models  include  the  Tresca  and  von  Mises  criteria,  which 
are  suitable  for  metal  plasticity,  and  the  Coulomb  and  Drucker-Prager  criteria, 
which  are  more  suitable  for  the  simulation  of  soil,  rock,  and  concrete. 

The  yield  stresses  in  both  the  Tresca  and  von  Mises  criteria  are  indepen¬ 
dent  of  pressure,  which  makes  these  models  unsuitable  for  simulating  the  pres¬ 
sure  dependent  material  behavior  exhibited  by  soil,  rock  and  concrete.  In 
contrast,  yield  stresses  in  the  Coulomb  and  Drucker-Prager  criteria  are  pres¬ 
sure  dependent.  For  more  information  on  these  models,  the  reader  should 
refer  to  Chapter  7  in  Owen  and  Hinton  (1980). 

An  elastic-plastic  strain-hardening  c^  model  was  implemented  into  JAM  to 
calculate  the  time-independent  skeletal  response  of  the  porous  solids.  The  cap 
model  enables  the  FE  code  to  model  nonlinear  irreversible  stress-strain  behav¬ 
ior  tmd  shear-induced  volume  changes.  Chapter  3  conuins  extensive  docu- 
menution  on  the  c^i  model. 


Element  Implemented  Into  JAM 

A  new  element  was  implemented  into  JAM  to  calculate  both  the  displace¬ 
ment  and  pore  fluid  pressures.  JAM  uses  an  eight-node  isoparametric  element 
with  16  displacement  and  four  pore  fluid  degrees  of  freedom.  Similar  eie- 
menu  were  used  by  Lewis  and  Schrefler  (1987).  Simon  et  at.  (1986a;  1986b). 
Zienkiewicz  (l98Sa),  Zienkiewicz  ft  at.  (1980),  and  Zienkiewicz  and  Shiomi 
( 1984).  Four  Gauss  integration  points  are  utilized  in  each  element. 
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Summary 


In  this  chapter,  the  work  of  Biot  and  other  investigators  was  briefly  de¬ 
scribed,  followed  by  a  discussion  of  the  modified  Biot  equations  implemented 
by  Lewis  and  Schrefler  and  modifications  that  were  made  to  those  equations. 
In  addition,  «]uations  were  derived  for  the  residual  forces.  Finally,  the  five 
constitutive  models  available  in  JAM  and  the  element  implemented  into  JAM 
were  described. 
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3  The  Cap  Model  and  its 
Implementation 


Introduction 

The  model  falls  within  the  framework  of  the  classical  incremental 
theory  of  work-hat* 'ening  plasticity  for  materials  that  have  time-  and  temper- 
ature-indqtendem  properties"  (Chm  and  Baladi  198S).  When  modelling  geo¬ 
logic  maten^la  subjected  to  stresses  ranging  from  one  to  several  hundred 
megapascals,  the  cap  model  has  several  desirable  features.  Of  primary 
importance  is  its  ability  to  model  volumetric  hysteresis  through  the  use  of  a 
strain-hardening  yield  surface  or  cap. 

In  this  chapter,  the  essential  features  of  the  cap  model  are  reviewed,  and 
the  steps  required  to  implement  the  cap  model  into  the  finite  element  code 
JAM  are  summarized.  After  a  brief  evaluation  of  the  loading  function  and 
flow  rule,  the  incremental  elastic-plastic  stress-strain  relations  are  outlined.  In 
addition,  the  cap  model  inqilemented  into  JAM  is  described,  and  the  equations 
are  developed  for  the  elastic-plastic  constitutive  matrix  and  the  plastic 
hardening  modulus.  Finally,  the  nunierical  inqilementation  of  the  cap  model 
itself  is  described.  The  reader  should  note  that  an  engineering  mechuiics  sign 
convention  is  used  in  which  compression  is  negative. 


Background 

The  cap  model  has  been  used  by  researchers  in  the  ground  shock  commu¬ 
nity  for  iq)proximately  20  years  to  siimilate  the  responses  of  a  wide  variety  of 
geologic  materials.  It  is  "predicated  on  the  fact  that  the  volumetric  hysteresis 
exhibited  by  many  geologic  materials  can  also  be  described  by  a  plasticity 
model,  if  the  model  is  based  on  a  hardening  yield  surface  which  includes  con¬ 
ditions  of  hyditswic  stress"  (Sandler  et  al.  1976).  The  model  was  first 
described  in  the  open  literature  by  DiMaggio  and  Sandler  (1971).  The 
FORTRAN  source  code  for  the  model  was  published  by  Sandler  and  Rubin 
(1979) 
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As  stated  above,  the  cap  model  has  several  desirable  features.  Of  primary 
importance  is  its  ability  to  model  volumetric  hysteresis  through  the  use  of  a 
strain-hardening  yield  surface  or  cap.  The  cap  model  may  also  be  formulated 
with  a  nonlinear  failure  surface,  with  linear  or  nonlinear  elastic  moduli,  or  as 
a  hinction  of  the  third  stress  invariant.  With  the  appropriate  selection  of 
material  parameters,  it  can  be  used  as  a  linear  elastic  or  linear  elastic-perfectly 
plastic  material  model.  Sandler  and  Rubin  (1979)  demonstrated  notable  fore¬ 
sight  with  their  use  of  function  statements  within  the  model,  which  allow  sub¬ 
stantial  changes  to  be  made  to  the  cap  model’s  potential  functions  with  little 
programming  effort. 

Modifications  and  expansions  of  the  original  model  were  made  by  several 
researchers.  Effective-stress  versions  of  the  cap  model  are  reported  in 
Baladi  (1979),  Baladi  and  Akers  (1981),  and  B^adi  and  Rohani  (1977,  1978, 
1979).  A  transverse-isotropic  cap  was  developed  by  Baladi  (1978)  and  an 
elastic-viscoplastic  cap  model  by  Baladi  and  Rohani  (1982).  Rubin  and 
Sandler  (1977)  developed  a  high-pressure  cap  model  for  ground  shock  calcu¬ 
lations  due  to  subsurface  explosions.  Baladi  (1986)  developed  a  "complex" 
strain-dependent  cap  model,  which  required  39  model  parameters,  for  ground 
shock  calculations  of  a  dry  cemented  sand.  In  addition,  several  versions  of 
the  cap  model  are  described  in  the  text  by  Chen  and  Baladi  (198S). 

In  formulating  the  cap  model,  DiMaggio  and  Sandler  (1971)  complied  with 
the  constraints  imposed  by  Drucker’s  stability  posmlate.  Drucker’s  stability 
postulate  is  sufficient,  although  not  necessary,  to  satisfy  all  thermodynamic 
and  continuity  requirements  for  continuum  models  (Sandler  et  al.  1976). 
Satisfying  Dmcker’s  stability  postulate  insures  uniqueness,  continuity,  and 
stability  of  a  solution  and  provides  a  mathematical  problem  that  is  properly 
posed.  Rubin  and  Sandler  (1977)  state  that  "...the  numerical  solution  to  a 
properly  posed  problem  can  proceed  without  the  fear  that  the  results  will  be 
strongly  dependent  on  errors  of  approximation  of  initial  and  boundary  condi¬ 
tions.  round  off  error,  etc."  Drucker  (19S1)  defines  a  work-hardening  mate¬ 
rial  as  one  that  remains  in  equilibrium  under  an  added  set  of  stresses  applied 
by  an  external  agency.  It  also  means  that  "(a)  positive  work  is  done  by  the 
extemr.I  agency  during  the  application  of  the  added  set  of  stresses  and  (b)  the 
net  work  performed  by  the  external  agency  ever  the  cycle  of  application  and 
removal  is  positive  if  plastic  deformation  has  occurred  in  the  cycle" 

(Drucker  1951)  Drucker  ( 1 950)  states  these  two  conditions  in  a  mathematical 
formal  as 

da^Jdt^J  >0  and  do,jdt^j  2:0 

The  first  sutement  constrains  a  model  such  thai  strain-softening  may  not 
occur.  The  second  sutement  unplies  (a)  the  loading  function  or  yield  surface 
must  be  convex  and  (b)  the  plastic  strain  increment  vector  must  be  normal  to 
the  yield  surface,  which  means  that  an  associated  flow  rule  must  be  used 
Ihese  are  the  constraiius  unposed  by  Drucker’s  subility  postulate. 
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Loading  Functions  and  Flow  Rule 


Drucker’s  criteria  for  stability  permits  considerable  flexibility  in  the  func¬ 
tional  forms  of  the  loading  function  /.  Since  Drucker’s  stability  postulate 
requires  the  yield  surface  and  plastic  potential  surface  to  coincide,  the  loading 
function  /  implicitly  represents  both  the  yield  and  potential  surfaces.  For  a 
perfectly  plastic  material,  a  general  form  of  the  losing  function  may  be  writ¬ 
ten  as 


f(Oij)  =0 


and  as 


/(ay-./c)=0  3.2 

for  a  strain-  or  work-hardening  material,  where  x  is  a  hardening  parameter 
that  acts  as  an  "..internal  sute  variable  that  measures  hardening  as  a  function 
of  the  history  of  plastic  volumetric  strain"  (Sandler  and  Rubin  1979).  For 
isotropic  materials,  the  loading  function  may  be  expressed  in  terms  of  stress 
invariants,  e.g.. 


where  7,  =  =  the  trace  of  the  stress  tensor  and 

Jjp  =  'As  jjSjj  -  the  second  invariant  of  the  deviatoric  stress  tensor. 
This  is  the  form  of  the  loading  function  used  in  most  versions  of  the  cap 
model.  The  loading  function  is  assumed  to  be  isotropic  and  is  comprised  of 
two  surfaces,  an  ultimate  failure  envelope  and  a  strain-hardening  surface  or 
cap.  The  failure  envelope,  which  is  fixed  in  space  and  symmetric  about  the 
hydrostatic  axis,  limits  the  maximum  shear  stresses  in  the  material  and  is  ex¬ 
pressed  as 


The  cap.  which  moves  as  plastic  deformations  occur,  is  represented  as 
/  -  H  ( ,  ^Jip  .  *  )  “  \^'^2D  ■  Ft  7|  ,*  )  =  0 

The  hardening  parameter  x  is  generally  taken  to  be  a  function  of  the  plastic 
volume  strain  (Chen  and  Baladi  198S) 
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*  = 


3.6 


I 

Equation  3.6  allows  the  cap  to  expand  and  contract.  By  allowing  the  cap  to  ^ 

contract,  one  can  limit  the  amount  of  dilation  that  a  material  may  develop 

when  its  stress  path  moves  along  the  failure  envelope  h.  This  form  of  the 

hardening  parameter  is  typically  used  for  soil-like  materials  that  do  not  exhibit 

significant  dilation  during  failure.  For  rock-like  materials,  the  hardening 

parameter  may  be  written  as  ^ 


In  this  form,  the  cap  is  only  permitted  to  expand,  thus  allowing  a  material  to 
dilate  while  its  stress  path  moves  along  the  failure  envelope,  i.e.,  when 
At  *  0.  Both  Equations  3.6  and  3.7  produce  hysteresis  during  an  imposed 
hydrostatic  load-unload  cycle  (Baladi  and  Akers  1981). 

The  plastic  loading  criteria  for  the  loading  fuiiction  /  are  given  by 


>  0  loading 
0  neutral  loading 
<  0  unloading 


3.8 


(Baladi  and  Akers  1981).  These  criteria  iirqply  that  during  loading  from  a 
point  on  a  given  yield  surface  a  stress  increment  tensor  do^j  (when  viewed  as 
a  vector)  will  poim  outward  (Rohani  1977).  Plastic  strains  will  only  occur 
under  this  condition.  During  unloading,  the  stress  vector  points  inward,  and 
the  material  will  behave  elastically.  Neutral  loading  occurs  when  the  stress 
vector  is  tangem  to  the  yield  surface.  During  neutral  loading,  no  plastic 
strains  are  produced  in  the  case  of  a  work-hardening  material  (Rohani  1977). 
This  is  referred  to  as  the  ‘contiiniity  condition*,  and  its  satisfaction  leads  to 
the  coincidence  of  elastic  and  plastic  constitutive  equations  (Chen  and 
Baladi  198S). 


Drucker  ( 19S1)  has  shown  that  the  plastic  strain  increment  tensor  for  a 
work-hardening  material  may  be  wrinen  as 


0 


if  /  -  0  and  ^  dffy  >  0 


do., 


d/ 


3.9 


if  /  <  0.  or  /  *  0  and  S  0 


which  is  identical  to  the  expression  used  for  elastic-perfectly  plastic  materiais. 
The  term  dX  is  a  positive  factor  of  proportionality  that  is  nonzero  only  when 
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plastic  defonnations  occur  (Baladi  and  Akers  1981).  For  the  cap  model,  the 
loading  function  /  may  take  the  form  of  either  Equation  3.4  or  3.S. 

Derivation  of  Incremental  Elastic-Plastic  Stress- 
Strain  Relations 

The  basic  premise  in  the  formulation  of  the  cap  model  and  all  elastic-plas¬ 
tic  constitutive  models  is  that  certain  materials  are  capable  of  undergoing 
small  plastic  (permanent)  strains  as  well  as  snull  elastic  (recoverable)  strains 
during  each  loading  increment  (Baladi  and  Akers  1981).  This  may  be  ex¬ 
pressed  mathematically  as 


3.10 


where  dt  jj  conq^nents  of  the  total  strain  increment  tensor, 

dt*j  components  of  the  elastic  strain  increment  tensor,  and 
dt^^j  »  components  of  the  plastic  strain  increment  tensor. 

This  equation  simple  states  that  the  total  strain  increment  is  equal  to  the  sum 
of  the  elastic  and  plastic  strain  increments.  In  its  most  general  form,  the 
elastic  strain  increment  tensor  may  be  expressed  as 


where  )  «  the  material  response  funaion,  which  may  be  a 

function  of  stress.  For  isotropic  materials,  the  elastic  strain  increment  tensor 
may  be  expressed  as 


dJ. 
9  k 


3,12 


where  »  o  -  ( 7,  /  3 )  6  *  the  deviatoric  stress  tensor. 

»  the  Kronecker  delu,  and 

K  and  G  are  the  elastic  bulk  and  shear  moduli,  respectively. 

The  elastic  bulk  and  shear  moduli  tnay  be  constants  or  functions  of  stress  or 
strain  invariants,  e.g., 

K  =  KiJ1.J2p.Jyp)  3  13 

G  *  ^iJ\.Jip.Jyp) 


where  Jyp  «  '/»s  *  the  third  invariant  of  the  deviatoric  stress 

tensor. 
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Chen  and  Baladi  (1985)  discuss  the  thermodynamic  restrictions  to  the  pos¬ 
sible  forms  of  the  above  equations.  Permissible  functional  forms  of  K  and  G 
must  not  generate  energy  or  hysteresis  and  must  maintain  the  path-independent 
behavior  of  elastic  materials.  Thus,  the  bulk  and  shear  moduli  should  be 
limited  to  the  following  forms  (Chen  and  Baladi  1985) 

K-KU„e^j)  3.14 


Inclusion  of  the  plastic  strain  tensor  into  the  functional  forms  of  K  and  G  is 
permitted  since  plastic  strains  are  constant  during  periods  of  elastic  deforma¬ 
tion.  Under  these  restrictions,  the  hydrostatic  and  deviatoric  components  of 
the  elastic  strain  increment  tensor  (Equation  3.12)  may  be  written  as 


d( 


e 

kk 


dJ^ 

3Ar(y,.*J) 


3.15 


and 


2C(yjo.<S) 


3.16 


Combining  Equations  3.9  and  3.12,  the  total  strain  increment  tensor  can  be 
written  as 


df 


y 


dJ, 


dSy 

2G 


dX 


K. 

do,, 


3.17 


The  plastic  strain  increment  tensor  (Equuion  3.9)  may  also  be  expressed  in 
terms  of  the  hydrostatic  and  deviatoric  cmnponenis  of  strain.  Applying  the 
chain  rule  of  differentiation  to  the  right-hand  side  of  Equation  3.9  results  in 


df 

d/l  dOy 


a/ 


d/Z 


2D 


dp. 


da. 


2D 


3.18 


which  simplines  to 


iLs 

dy. 


a/ 


2^J2D  ^p2D 


'ij 


3.19 


Multiplying  both  sides  of  Equation  3.19  by  gives  an  expression  for  the 
plastic  volumetric  strain 
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dd,  =  3  d\ 


By  definition,  the  deviatoric  component  of  the  plastic  strain  increment  tensor 
is  written  as 

defj  -  3-2' 


Substitution  of  Equations  3.19  and  3.20  into  Equation  3.21  gives 

deP.^-J!h _ ^s,^d\lL 


2  p2D  ^^2 


d5,, 


The  proportionality  factor  d\  must  be  determined  prior  to  evaluating  any  of 
the  above  plastic  strains.  Baladi  and  Akers(1981).  Chen  and  Baladi(198S), 
and  Rohani(1977)  outline  the  methods  required  to  evaluate  the  proportionality 
factor.  Those  methods  are  included  here  for  conpleteness. 

Using  Equations  3.4,  3.S,  and  3.6  or  3.7,  the  toul  derivative  of  the 
loading  function  /  may  be  expressed  as 

df  *  *  - ! - tL —  s.jdSj:  *  ^  «  0  71 


2JJ20 


This  expression  is  known  as  the  "consistent  condition"  for  strain-hardening 
materials  (Chen  and  Baladi  1985).  Using  Equations  3.15,  3.16,  and  3.20, 
Equation  3.23  may  be  manipulated  to  give 


Cde*. 


ho 


a.  .3  24 
-  “  0 


Substitution  of  Equation  3. 10  into  Equation  3.24  produces 

3  Kidt/^^  -  )  -JY  *  ij  ~  )  — p== 


^  3,1, 
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The  stress  increment  tensor  may  be  written  as 

doij  =  Kd€t^6ij  +  IGde^j 


kk 


3/ 


p2D 


9K 


3/ 

2 

+  G 

3/ 

'.33/ 

3/ 

ax 

ay, 

dp2D 

3y, 

ax 

3eu 

6. 

3/ 

.  1^- 

ay, 

^p2D 

■ij 

3.29 


Equations  3.28  and  3.29  are  the  general  constitutive  equations  for  an  elas¬ 
tic  work-hardening  plastic  isotropic  material  (Chen  and  B^adi  1985).  To  use 
these  equations,  one  must  first  define  the  loading  function  /,  the  functional 
forms  of  the  elastic  moduli  ilT  and  G,  and  the  hardening  parameter  k  for  the 
material  of  interest. 


ElastiC'Plastic  Constitutive  Matrix 

In  the  following  section,  the  equations  for  the  elastic-plastic  constinjtive 
matrix  are  formulated.  The  equations  are  written  in  matrix  format  to  render  a 
more  compact  form  of  the  equations.  The  development  of  the  elastic-plastic 
constitutive  matrix  follows  the  derivuion  of  Owen  and  Hinton  (1980). 

The  loading  or  yield  function  /  for  a  general  work-  or  strain-hardening 
elastic-plastic  model  (Equation  3.2)  may  be  rewritten  as: 

/(*,«)  -  F{0)  -  *(«)  -  0  ^  ^ 

where  (in  matrix  format)  «  is  the  ve^or  of  normal  and  shear  stresses  and  x  is 
the  hardoiing  parameter  that  controls  the  expansion  of  the  yield  surface. 

Recall  that  in  the  cap  model,  the  cap  itself  is  the  only  strain-hardening  yield 
surface;  the  failure  envelope  is  not  a  hardening  surface.  Equation  3.^  may 
be  differentiated  to  give: 

df  •  —d»  *  ^dn  “  0  3.31 

d«  dx 


or  in  another  form: 
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a^do  -  A  d\  =0 


where 

da  dOii 


A  ^  -J-  2L  dK 
d\  dx 


Owen  and  Hinton  refer  to  the  vector  a  as  the  flow  vector.  The  scaler  A  will 
be  identified  as  the  plastic  hardening  modulus. 

The  total  strain  increment  tensor  (Equation  3. 10)  may  be  wrinen  in  matrix 
format  as 

dt  «  dt‘  *  dtP 

By  substituting  for  both  the  elastic  and  plastic  strain  increments,  i.e.,  using 
the  matrix  equivalents  of  Equations  3.9  and  3.11,  the  following  expression  is 
obtained: 

dt  -  D'^de  *  d\K  3.36 


where  D  is  the  matrix  of  elastic  material  constants  and  the  inverse  of  the 
material  response  function. 

After  multiplying  Equation  3.36  by  a^D  one  obtains; 

aJ Ddi  -  ♦  aJOad'k 


which  may  be  refined  further  by  eliminating  a^de  with  the  use  of  Equa¬ 
tion  3.32  to  produce: 

a^Ddt  -  I  A  *  a^Da  [dX 

This  leads  to  an  expression  for  the  scaler  term  dk: 

a^Ddt 

d\  •  - - 

A  ♦ a^Da 


This  teim  gives  the  mi^tude  of  the  plastic  strain  incremem  vector  and  is  the 
matrix  form  of  Equation  3.27.  Note  Uuu  in  the  text  by  Owen  and 
Hinton  (19M).  this  eiqireuion  was  printed  incorrectly. 
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Substituting  Equation  3.44  into  EquatitMi  3.34,  substituting  for  e^,  and  rear¬ 
ranging  produces: 

A  =  -KIJLM.  3.4 

3k  3<r 

The  plastic  hardening  modulus  A  will  be  dependent  upon  the  functional  form 
of  the  loading  or  yield  function  /  and  the  hardening  function  used  in  the  cap 
model. 


Cap  Model  Implemented  into  JAM 


The  following  section  describes  the  version  of  the  cap  model  implemented 
into  the  fmite  element  code  JAM.  The  functional  forms  of  the  equations  are 
outlined,  and  the  cap  model  is  described  in  more  detail. 


Two  elastic  response  functions  govern  the  behavior  of  the  model  in  the 
elastic  regime.  The  elastic  bulk  modulus  is  defined  by  the  following  equation 


3.46 


where  AT,  ^  the  initial  elastic  bulk  modulus  and 
AT,  and  Ar2  are  material  constants. 

The  elastic  bulk  modulus  prescribes  the  unloading  moduli  in  pressure-volume 
space.  The  three  material  constants  may  be  determirnd  from  the  unloading 
data  obuined  during  hydrostatic  loading  tests.  The  elastic  shear  modulus  is 
defined  by  the  following  equation 


G(J. 


2D 


*{*) 


Cj  exp(  G2  f^if) 
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where  G,  »  the  initial  elastic  shear  modulus  and 
G]  and  G2  are  material  constants. 

The  elastic  shear  modulus  prescribes  the  unloading  moduli  in  principal  stress 
difference-principal  strain  difference  space  The  three  material  constants  may 
be  determined  from  the  unloading  shear  data  of  triaxial  compression  tests. 

The  constaius  for  the  elastic  bulk  and  shear  moduli  may  also  be  determined 
from  uniaxial  strain  con^ression  tests,  i.e.,  from  the  unloading  slopes  of  the 
stress  path  ( *  2 G/ AT)  and  stress-strain  curves  ( -  AT  ♦  4/3  G). 


In  the  current  model,  the  failure  envelt^  portion  of  the  loading  function  / 
is  defined  by  a  modified  Orucker-Prager  failure  surface  (see  Figure  3. 1)  of  the 
form 
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Pkflm 


h  (Ji, 


RDiptk  Cap 


Cap  _ 

TraoaitMO 


»«)  •  U«) 


3.1.  Cap  model  yield  surfaces 


h{Ji  ,p2D  )  ”  p2D  ~  ~  Cexp(By,  )|  if  y,  >!,(«) 

where  ,  S  and  C  are  material  consunis.  These  constants  may  be  determined 
ftom  the  locus  of  triaxial  compression  failure  dau  ploned  in  Uk  ai^ropriate 
stress  space.  The  strain-hardening  yield  surface  or  cap  is  described  by  the 
following 

H(Ji 

-  I([X(IC)  -  L(ic)]^  -  [y,  -  ify,  <L{k) 


where  X(  k  )  and  L(k)  define  the  values  of  Jf  at  the  intersection  of  the  cap 
with  the  Ji  axis  and  at  the  center  of  the  cap.  respectively  (see  Figure  3. 1 ); 

K  is  the  hardening  parameter,  which  is  equal  to  the  plastic  volumetric  strain. 
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Equations  3.48  and  3.49  also  indicate  the  value  of  Jj  determines  which  of  the 

two  yield  surfaces  should  be  used.  /?  is  the  ratio  of  the  major  to  minor  axes 
of  the  elliptical  cap  and  has  the  following  hinctional  form 

R[L(k)]  =  /?,.  >  /?i[  1  -  exp(/?2{^(x)  -  I-o})J 


where  /?2  and  Lq  are  material  constants;  Lq  defines  the  initial  location 

of  the  cap. 

Chen  and  Baladi  (1985)  explain  that  the  functional  fonn  of  the  cap  was 
chosen  such  that  the  tangent  of  its  intersection  with  the  failure  envelope  is 
horizontal.  This  condition  is  guaranteed  by  the  following  relationship  between 
X(k)  and  L(  k  ) 

X(k)  -  L(k)  -  Rhil(K),^)  3.52 


where 


1  0if/(it)^0 


The  hardening  function  for  tii's  model  is  defined  by 


which  may  be  rewritten  in  the  following  form 


1  ^  kk 

X(  K )  «  -L  In  _  ♦  1  ♦  Xo 

D  W  ® 


where  D.  W  and  Xq  are  material  constants.  W  esublishes  the  maximum 
plastic  voliunetric  strain  the  material  can  develop;  Xq  .  like  Lq  .  defines  the 
initial  location  of  the  cap. 

As  described  previously.  Orucker’s  stability  posnilaie  places  limits  on  the 
functional  forms  of  the  equations  in  the  cap  model.  Sandler  and  Rubin  ( 1979) 
specify  some  of  those  limitations  ,  (a)  Q{  )  must  decrease  monotonically 
with  increasing  values  of  ;  (b)  to  avoid  work-softening,  the  functions 
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X(k)  and  L(  k  )  must  be  continuous  and  monotonicaily  increasing  functions 
and 

4^2:0  and  |^<0 
dJi  dx 

(c)  the  cap  must  extend  from  the  Jj  axis  to  a  point  on  or  below  the  failure 
envelope  h,  i.e., 

F[X(k),k]^0 


and 


(d)  within  the  yield  surfaces  defined  by 

<  QiJO  fory,  >L(k) 


3.57 


3.58 


and 


p2D  for  Hk)  ^J^ttXiK) 

the  material  response  must  be  isotropic  elastic.  Sandler  and  Rubin  (1979) 
explain  that  if  the  itwquality  in  Equation  3.57  is  true,  then  a  gap  exists  be¬ 
tween  the  cap  f/  and  the  failure  envelope  h  and  a  von  Mises  type  failure  sur¬ 
face  is  used  as  a  transition  between  the  two  yield  surfaces  (Figure  3.1).  The 
yield  surface  for  itL{  k  )  is  thus  defined  by  the  following  expression 

-min(f(L(«).*l  .  C(y,)) 


In  the  evaluation  of  plastic  hardening  modulus  A.  one  need  only  be  con¬ 
cerned  with  the  functional  form  of  the  cap  since  it  is  the  only  hardening  sur¬ 
face  in  the  model.  In  order  to  evaluate  A,  each  of  the  three  terms  in  Equation 
3.45  must  be  evaluated.  Recalling  Equation  3.50.  the  first  term  may  be  ex¬ 
panded  as 

dH  ^  dHdX  ^  dH  dX 

a«  ‘  ■Jxd*  “  ^  ^ 
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Each  of  these  terms  is  evaluated  as 


Z)U{* .  IV 


=  2(1  -X) 


Combining  the  terms  gives 


aw  _  2(1-20 

ok  ♦ 


The  second  term  in  the  expression  for  A  is  evaluated  as 


since  k  «  .  The  fuial  term  in  ^4  can  be  expanded  as 


aw  aw.  ^ij  aw 
da,,  ay,  -/7 — ./y— 


Recalling  that  >  0  and  substituting  Equations  3.64-3.66  into  Equa¬ 
tion  3.4S,  one  gets 

.  6(X  -  L)  dH 

A  B  - - - 3 

*  W)^\ 


as  an  expression  for  A.  Simplifying  this  further  by  evaluating 

42 . 2u,  -  i) 


and  substituting  into  Equation  3.67.  one  obuins  the  final  expression  for  the 
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plastic  hardening  modulus 


^  \2{X-L)(Jy-L) 

D(<U  •  O') 


Numerical  Implementation  of  the  Cap  Model 

introduction 

The  numerical  algorithm  for  the  cap  model  was  published  by  Sandier  and 
Rubin  in  an  attempt  *to  facilitate  the  general  use  of  the  cj^  model  in  dynamic 
compuutions,  as  well  as  in  model  fitting"  (Sandier  and  Rubin  1979).  The  cap 
model  algorithm  was  designed  for  use  in  either  finite  element  or  finite  differ¬ 
ence  codes  and  is  a^iicabie  to  both  static  and  dynamic  problems  (Chen  and 
Baladi  1985).  Of  notable  foresight  on  the  part  of  the  designers  was  their  use 
of  fimction  suteirxnts  within  the  model,  which  allow  substantial  changes  to  be 
made  to  the  cap  model’s  potential  functions  with  little  programming  effort. 
This  feature  has  allowed  investigators  to  simulate  a  wide  variety  of  natural  and 
man-made  materials  with  high  degrees  of  fidelity  between  model  and  material 
response.  Despite  the  many  published  variations  of  the  cap  model,  the 
original  cap  model  algorithm  developed  by  Sandler  and  Rubin  still  forms  the 
foundation  of  most  curreiu  cap  model  algorithms. 

The  cap  model  algorithm  is  essentially  an  implementation  of  Equa¬ 
tion  3.29.  To  march  the  calculation  through  time,  the  user  must  input  the 

stresses  and  the  location  of  the  cap  at  time  /,  which  is  explicitly  defined 

by  the  term  1(’k)  xni  implicitly  defined  by  the  hardening  parameter  ,  and 
the  strain  increments  from  the  solution  of  the  field  equations  for  the  current 

time  step  The  cap  model  returns  the  new  stresses  and  the 

updated  cap  location  and  hardening  parameter  /(  )  and  at  time 

t*  At.  A  given  strain  incremeM  may  invoke  four  different  types  of  stress 
paths  that  coincide  with  four  different  algorithms  within  the  cap  model  itself; 

(a)  an  elastic  algorithm. 

(b)  a  failure  envelope  algorithm. 

(c)  a  hardening  cap  algorithm,  or 

(d)  a  tension  cutoff  algorithm. 

In  the  following  text,  a  description  of  the  four  numerical  algorithms  is  provid¬ 
ed  The  descriptions  are  basrj  upon  previous  descriptions  by  Baladi  and 
Akers  (1981).  Q»en  and  B»'adi  (1985).  Sandler  vid  Rubin  (1979),  and  Meier 
(1989).  To  simplify  the  p  natation,  a  description  of  the  cap  model's  re¬ 
sponse  in  the  tensile  regime  will  be  deferred  to  the  later  pan  of  this  section. 
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Elastic  algorithm 


To  start  the  numerical  procedure,  it  is  assumed  that  the  given  strain  incre¬ 
ments  produce  an  entirely  elastic  stress  path.  A  set  of  elastic  trial  stresses  are 
calculated  h-om 


%  =  %  ^3K  3.7( 

and 

.  'jy  .  2G  3.71 

The  elastic  trial  stresses  are  tested  with  respect  to  the  tension  cutoff,  the  fail¬ 
ure  envelope,  and  then  the  c^.  If  these  surfaces  are  not  violated  by  the  trial 
stresses,  the  actual  stress  path  is  an  elastic  path,  and  the  new  stresses  are  the 

elastic  trial  stresses,  i.e.,  -  ^7,  and 


Failur*  anvalopa  algorithm 

If  the  following  conditions  exist  when  the  elastic  trial  stresses  are  tested 
with  respect  to  the  failure  envelope, 

•  /  ^■^2D  )  *  /  ■  C  ( ^  ® 

then  the  elastic  trial  stresses  have  violated  the  failure  envelope,  and  the  given 
strain  increment  must  be  a  combination  of  elastic  and  plastic  strains.  The  trial 
stresses  must  be  corrected  such  that  (a)  tlw  fmal  stress  state  falls  on  the  failure 
surface  and  satisfies  the  following  relation 

A('*^'7,  ,^'*^'720  )-0 

and  (b)  the  resulting  elastic  and  plastic  strain  increments  add  up  to  the  given 
strain  increments  ,y . 

The  mathematical  statement  that  requires  the  fmal  stresses  to  lie  on  the 
fixed  failure  surface  is  given  as 

dh  -  4^dai.  -  0  3.73 

Assuming  small  strain  increments.  Equation  3.73  can  be  numerically  i^roxi- 
mated  by  the  following  expression 
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d/j  =  J 


'ho  -  Qi^Jx)  ^ 


which  reduces  to 


dh  =  t/%n  -  Q(  %) 


since  the  fmal  stress  point  must  lie  on  the  failure  surface,  i.e., 

-  Q(  -  0  - 

Equation  3.7S  may  be  substituted  into  Equation  3.73  and  expanded  in  the 
following  manner 


dh  dh 


dJL 


UB_  dOii  3.77 
don 


i^dy,  *  -j _ 

‘  8/^ 


where  dy]  *  ^yj  -  Vj  and  dsjj  •  -  'j,y.  From  Equations  3.70  and 

3.71,  we  know  that  dJ^  «  and  ds^j  «  2G'*^'de,y,  and  these 

expressions  may  be  substituted  imo  Equation  3.77  to  give 

-e(^y,)- 3* 

3.78 

♦  ^  '*^'de 

I -  F= 

V-^20  ^  V2D 
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An  expression  for  may  be  developed  in  the  following  manner 

'•"y,  .V,  OAT*;, 

■  v,  3*1 

.  ^y,  -  3Ar*f, 

where  is  defined  by  Equation  3.80.  A  ’temative*  value  of  may 

be  calculated  from  Equation  3.81;  this  value  is  tentative  because  it  must  be 
tested  against  the  current  position  of  the  cap.  which  is  defined  by  the  value  of 

If  '*^'7]  <  £.('«),  which  indicates  the  stress  poim  has  violated  the  cap, 
then  comer  coding  is  required,  i.e..  the  cap  must  intersect  the  failure  envelope 

forming  a  comer,  and  the  value  of  '*^'71  tmisi  be  adjusted.  Adhering  to  the 
imposed  conditions  of  notmality.  a  stress  state  lying  on  the  failure  envelope 
pr^uces  dilatant  plutic  volumetric  strains.  Since  cap  expansion  can  only 
result  from  compressive  plastic  volumetric  strains,  the  cap  is  stationary,  and 
the  new  stress  state  can  not  move  beyond  the  interseaion  of  the  cap  and  the 

failure  envelope.  Thus,  the  final  stress  state  is  '*^'7]  -  L(  'k  ).  and  the  up¬ 
dated  hardening  parameter  is  '*^'k  -  'k. 


Ch«pi*r  3  The  Cap  Modal  and  Itt  Implamantaiion 


If  >  L{^k),  then  the  final  stresses  will  depend  upon  the  form  of 

the  hardening  function.  Equation  3.7  is  the  simplest  form  of  the  hardening 
function  to  use  because  it  only  permits  plastic  volumetric  compaction,  i.e.,  the 
cap  is  only  allowed  to  expand.  As  in  the  above  case,  a  stress  state  lying  on 
the  failure  envelope  produces  dilatant  plastic  volumetric  strains.  Since  the 
hardening  function  defined  by  Equation  3.7  prescribes  no  cap  movement  due 
to  dilatant  volumetric  strains,  the  cap  is  stationary.  Thus,  the  final  stress  state 

is  ,  i.e.,  no  adjustment  is  required,  and  the  updated  hardening  parame¬ 
ter  is  -  'k.  If  the  hardening  function  takes  the  form  of  Equation  3.6, 

which  allows  the  c^  to  expand  and  contract,  the  cap  is  adjusted  (in  this  case 
contract^)  to  a  position  prescribed  by 


f*Af/  _ 


'/ 


3.82 


and  a  tentative  value  of  is  obtained.  The  new  position  of  the  cap  must 
be  compared  to  the  value  of  .  If  the  cap  has  contracted  such  that 

'*^'7j  <  ==  both  and  '*^'7,  must  be  adjusted 

such  that  This  is  accomplished  by  starting 

with  the  following  relation 


^J^  -  3Kdel^  *  -  ‘I  ♦  — 


3.83 


'/ 


eliminating  dth^,  which  is  the  third  unknown,  by  substituting  the  following 


d^kk 


^J,  -  '/ 


dl 

*  3K 

^*kk 

'1 

3.84 


from  which  one  can  show  that 
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i 


•  • 


t 


i 


t 


1 


'"^7,  =  -  3a: 


-  7 


dl 

*  3K 

^4k 

•'1 

dl 

+  3a: 

-  3K(^Ji  -  7) 

dl 

+  3a: 

>1 

which  in  mm  sin^lifies  to 


t*Ati  . 


n 


%  *  3K  7 


.1 


*  3K 


Having  calculated  the  final  value  of  ,  we  must  calculate  the  new  com¬ 
ponents  of  the  deviatoric  stress  tensor  The  expressions  for 

are  develqied  below. 

Recognizing  the  path  independence  of  linear  elastic  constimtive  equations, 
we  can  write 


Substimting  Equation  3.71  imo  Equation  3.87  and  performing  a  simple  manip¬ 
ulation  one  obtains 

'*^7  .  -  2Gde^ 

*  ij  *ij 


Recalling  Equation  3.22  and  recognizing  that  — 1—  -  1 .  we  can  write 

& 


■ij  “ 
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which  can  be  substituted  into  Equation  3.88  giving 


=  ^s-  - 


After  rearranging  the  above  equation  one  obuins 


Squaring  each  side  of  Equation  3.91  and  multiplying  by  'A  produces 


l*Ali  I  . 

hD  *  * 


Taking  the  square  root  of  each  side  and  rearranging  terms,  one  obtains 


Replacing  the  right-hand  side  of  Equation  3.93  with  the  expressions  in  Equa¬ 
tion  3.91,  one  obtains 


which  may  be  rewritten  for  our  use  as 


no  E, 


to  calculate  the  new  deviator  stress  tensor  components. 
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Cap  algorithm 


If  the  failure  envelope  is  not  violated  by  the  elastic  trial  stresses,  the  trial 
stresses  are  checked  against  the  loading  function  for  the  cap.  If  the  following 
conditions  exist 


V,  <X(‘k) 


then  the  c^  algorithm  is  invoked  and  the  position  of  the  cap  is  adjusted  until 


.1 


‘k)  =  0 


An  iterative  procedure  is  used  in  the  cap  algorithm.  To  start  the  proce¬ 
dure,  a  trial  value  of  dl^‘^  is  assumed  in  order  to  calculate  a  new  trial  cap 
position  =  /<')  *  '/  +  where  the  superscript  /  denotes  an 

iterative  value.  In  addition,  trial  values  of  x  ^ ,  L  ( x  ^ ,  X ( x  ^ .  and 

are  computed.  Finally,  a  trial  value  of  Jj  is  computed  from  the 
following  relation 

3  ^ 

If  £1  X(  x^'^ ),  a  smaller  value  of  dl^'^  is  assumed.  If 

y/*'  a  L(  x^'^ ),  a  larger  value  of  dl^'^  is  assumed.  This  process  is  carried 

on  until  the  condition  L(  x^'^ )  sy/'^  <  X(  x^'^ )  is  satisfied.  The  final  value 
of  /  is  one  which  satisfies  the  following  equation  to  some  desired  accuracy 

.fTTir.  .  _  .rE7~  3.97 


where 


^•^1  /<•>  ,(•) 
Jt  » « 


The  derivation  of  Equation  3.97  is  outlined  in  the  following  text. 
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If  we  start  with  Equation  3.22  and  substitute  for  d\  using  Equation  3.20, 
one  obtains 


1  dH  ],dH 


(ipiD 


Recognizing  that  f  ^  H,  d///3^72D  ~  ^  dHIdJ^ 
we  can  rewrite  Equation  3.99  as 


=  -dF/dJ^  = 


defj  . 


3.100 


Substituting  the  above  into  Equation  3.88  and  rearranging  gives 


c*{, 


3  ^  p2D 


Performing  the  same  operations  on  the  above  equation  as  was  used  on  Equa¬ 
tions  3.91-3.93,  one  obtains  Equation  3.97. 

The  solution  of  Equation  3.97  is  obtained  through  the  use  of  an  iterative 
convergence  routine  known  as  the  modified  reguU  falsi  method  (Sandler  and 
Rubin  1979).  A  dimensionless  function  P{1)  is  defuied  as 


/(«)  ~jr 

/(«)  -  X(i() 


if 


I  ♦Aly  (0 
•'2D 


Gdel, 


El  .  Jl*AljU)  ^ 
ho  ^  V  •'2D  ^  —JJ— 


if  X((()<-^,  </.(«) 


x(«)  -y,'" 

L(k)  -  X(k) 


if  y/"iL(K) 


3.102 


where  the  solution  Pit)  -  0  is  also  the  solution  of  Equation  3.97.  If  we  can 
show  that  ^y,  <  ‘*^‘1  <  '/.  then  P(l)  is  bounded  and  monotonic  in  the 

strict  sense,  and  the  solution  P(  / )  «  0  is  unique  and  can  be  found  to  any 
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desired  degree  of  accuracy  (Sandler  and  Rubin  1979).  An  expression  for  the 
degree  of  accuracy  or  tolerated  enor  is  given  by 


<  NQ[X{k)] 


3.103 


where  a  tolerance  of  N/h{oo  ,  )  =  10  ^  (in  dimensionless  format)  is 

typically  used. 


To  show  that  <  ‘*^‘1  we  must  recognize  that  <  L  (  V ) 

must  be  true,  since  it  is  a  condition  for  invoking  the  cap  algorithm.  In  addi¬ 
tion,  since 


3.104 


we  know  that  ,  because  the  plastic  volumetric  strain  increment  is 

negative  during  volumetric  conduction,  i.e.,  when  the  ct^)  expands.  This 

means  that  the  fmal  value  of  must  lie  in  the  range 

<L('«) 

Now  let  us  determine  the  lower  limit  of  which  will  lead  us  to  the 

lower  limit  of  The  cap  exhibits  its  furthest  expansion  when  ^7j  is 

at  the  intersection  of  the  oq)  and  the  failure  envelope.  When  the  cap  is  in  this 

position,  »  ^y| .  which  implies  that  the  lower  bound  of 

is 

t*sti  .  »  ^y, 

The  upper  bound  of  is  simply  the  value  at  time  f.  i.e.,  '*^'1  »  '/. 

Combining  these  expressions,  the  range  of  must  be 

^y,  <  <  '/ 

With  the  above  conditions  satisfied,  the  solution  to  Equation  3.97  may  be 
obtained.  This  concludes  the  description  of  the  cap  algorithm.  A  description 
of  the  cap  model's  response  in  the  tensile  regime  follows. 

T«n«ii«  algorithm 

Sandler  and  Rubui  (1979)  recognized  that  soil  tensile  dau  is  seldom  ob¬ 
tained  in  the  laboratory  and  therefore  dealt  with  lensile  behavior  in  a  simple 
manner.  They  also  cautioned  potential  users  of  the  siiTd>listic  n^re  of  tlu 
cap  model  in  the  tension  regime.  A  tension  failure  response  is  invoked  if 
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^7)  >  r,  where  T  is  the  tension  cutoff  or  limit.  Sandler  and  Rubin  recom¬ 
mended  that  the  final  stresses  be  defmed  as  =  T  and  =  0 

when  the  tension  cutoff  is  exceeded.  For  materials  using  the  definition  of  the 
hardening  parameter  defined  by  Equation  3.6,  the  plastic  volumetric  strain  is 
defined  by 


=  d( 


kk 


%  -  % 


3.105 


and  an  updated  hardening  parameter  is  determined.  If  the  tension  cutoff  is  not 

exceeded  but  ^7]  >  0,  i.e.,  the  elastic  trial  stress  still  lies  in  the  tension 
regime,  the  stress  state  must  be  checked  against  both  the  failure  envelope  and 
the  von  Mises  transition  using  the  following  inequality 

S  min(G(^7,),  f[L(  'x),'*]) 


Stress  states  violating  the  von  Mises  transition  must  be  remmed  to  that  surface 
using  the  same  logic  implemented  for  the  failure  envelope.  Stress  states  lying 
on  the  von  Mises  transition  will  produce  no  plastic  volumetric  strains  due  to 
the  inqjosed  normality  conditions.  In  addition,  the  von  Mises  transition  is 
fixed  because  the  cap  hardening  surface  does  not  expand. 


Implementation  of  Cap  Model  into  JAM 

Two  basic  operations  that  are  associated  with  elastic-plastic  material  mod¬ 
els  must  be  performed  in  most  inqilicit  finite  element  codes;  (1)  the  construc¬ 
tion  of  the  elastic-plastic  constitutive  matrix  and  (2)  the  calculation  of  the 
residual  forces.  The  purpose  of  this  section  is  to  explain  how  these  two  oper¬ 
ations  were  affected  by  the  implementation  of  the  cap  model.  The  later  opera¬ 
tion  will  be  considered  first  since  it  is  a  straight  forward  process. 

After  the  strain  increments  at  time  Ar  are  obuined  from  the  solution  of 
the  field  equations,  the  stresses  at  time  t*At  in  each  element  are  calculated. 
The  residual  forces  at  time  r  >  Ar  are  then  calculated  based  on  the  stress  states 
in  the  elements.  In  this  operation,  no  substantial  changes  are  required  in  the 
cap  model  subroutines;  hence,  the  implementation  is  rather  simple 

However,  three  components,  the  elastic  constinitive  matrix  D,  the  plastic 
hardening  modulus  A ,  and  the  flow  vector  a,  are  required  to  calculate  the 

elastic-plastic  constitutive  matrix  The  calculation  of  the  elastic  constitu¬ 
tive  matrix  is  sinqile  and  needs  no  further  discussion.  In  JAM,  a  modified 
cap-model  subimitine  calculates  and  returns  values  for  the  flow  vector  and  the 
plastic  hardening  modulus.  This  subroutine  first  determines  which  of  four 
possible  regions  or  surfaces  a  stress  point  resides  in  or  on,  i.e.,  an  elastic 
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region,  the  failure  surface,  the  cap,  or  the  tension  cutoff.  The  plastic  harden¬ 
ing  modulus  is  nonzero  only  when  the  stress  point  falls  on  the  cap  since  it  is 
the  only  hardening  surface.  In  this  case,  the  plastic  hardening  modulus  is 
calculated  using  Equation  3.69.  The  flow  vector  is  calculated  by  numerically 
evaluating  Equation  3.66  when  the  stress  state  lies  on  the  failure  surface,  the 
cap,  or  the  tension  cutoff. 

To  complete  the  implementation,  one  must  provide  the  model  access  to  the 
material  constants  and  an  array  to  store  the  location  of  the  cap  for  each  nu¬ 
merical  integration  point. 


Verification 

To  insure  that  the  cap  model  was  correctly  incorporated  into  the  finite 
element  program  JAM,  several  laboratory  stress-  and  strain-path  tests  were 
numerically  simulated.  These  calculations  were  compared  to  the  output  from 
a  cap  model  driver  (Chen  and  Baladi  198S)  exercised  over  the  same  laboratory 
stress  and  strain  paths.  The  two  programs  should  produce  similar  if  not 
identical  results. 

The  following  tests  and  strain  paths  were  simulated:  a  hydrostatic  compres¬ 
sion  test  with  one  load/unload  cycle,  a  set  of  constant  radial  stress  triaxial 
compression  tests,  a  set  of  constant  mean  normal  stress  tests,  a  uniaxial  strain 
(Kg)  test  with  one  load/unload  cycle,  and  finally  a  test  with  a  K,,  load/constant 
axial  strain  (BX)  unload  cycle.  To  simulate  the  tests  with  the  finite  element 
program,  a  single  element  was  loaded  under  the  appropriate  boundary  condi¬ 
tions.  Several  loading  increments  were  utilized  during  each  calculation,  and  a 
convergence  tolerance  of  1  percent  was  satisfied  at  the  end  of  each  increment. 
The  output  at  the  end  of  each  increment  is  represented  by  a  symbol  on  the 
comparison  plots. 

The  simulated  hydrostatic  loading  test  consisted  of  an  applied  loading  to  a 
pressure  level  of  250  MPa,  followed  by  an  unloading  to  zero  pressure.  This 
calculation  exercised  the  logic  and  code  affecting  both  the  cap  movements  and 
the  elastic  algorithms  within  the  model  and  finite  element  program.  The  finite 
element  and  cap  model  driver  results  are  compared  in  Figure  3.2.  Output 
from  the  finite  element  program  matches  the  cap  model  driver  with  no  notice¬ 
able  errors. 


Four  constant  radial  stress  triaxial  compression  tests  at  radial  stresses  of 
25,  50,  1(X),  and  150  MPa  were  simulated.  Loading  was  terminated  prior  to 
reaching  the  ultimate  failure  surface,  anJ  unloading  results  were  acquired  for 
only  three  of  the  four  calculations.  The  values  of  principal  stress  difference 
calculated  by  the  finite  element  program  were  less  than  those  of  the  cap  model 
driver  (Figure  3.3).  The  magnitude  of  the  errors  decreased  when  a  larger 
number  of  increments  was  applied. 
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Figure  3.2.  Simulated  hydrostatic  loading  stress-strain  response 


vf  I  S 


Figure  3.3.  Simulated  triaxial  compression  stress-strain  response 

Four  coiatam  mean  noraiai  stress  tesu  at  confuung  pressures  of  23.  SO, 
100,  and  ISO  MPa  were  simulated.  As  with  the  simulated  triaxial  coropres- 
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Figure  3.4.  Simulated  stress-strain  response  of  constant  mean  normal 
stress  tests 

sion  tests,  loading  was  terminated  prior  to  reaching  the  ultimate  failure  sur¬ 
face;  unloading  results  were  not  acquired.  The  finite  element  and  cap  iiKXlel 
driver  results  are  conq)ared  in  Figure  3.4.  Errors  in  the  stresses  calculated  by 
the  finite  element  program  were  less  than  those  in  the  simulated  triaxial  com¬ 
pression  tests  (Figure  3.3). 

The  simulated  test  consisted  of  an  applied  loading  to  a  vertical  strain 
level  of  20  percent,  followed  by  an  unloading  to  a  small  value  of  vertical 
stress.  In  this  calculation,  the  finite  element  simulation  was  conducted  with  a 
displacement  controlled  boundary  condition.  This  type  of  loading  should 
produce  an  exact  match  between  the  finite  element  program  and  the  cap  model 
driver  since  no  strain  increment  iterations  are  required  in  the  finite  element 
program.  The  calculated  Kg  stress-strain  response  is  ploned  in  Figure  3.5  and 
the  stress  path  response  in  Figure  3.6.  There  are  no  noticeable  differences  be¬ 
tween  the  two  calculations. 

A  Ko  load/BX  unload  test  was  also  simulated  with  a  displacement  control¬ 
led  boundary  condition.  The  test  consi^  of  an  applied  loading  to  20  percent 
axial  strain,  followed  by  an  unloading  to  a  small  value  of  radial  stress  (Fig¬ 
ure  3.7).  In  this  calculation,  the  comer  coding  of  the  cap  model  was  exer¬ 
cised  as  the  stress  path  unloaded  along  the  failure  envelope  (Figure  3.8).  The 
calculated  resulu  suggest  a  proper  implementation  of  this  logic. 
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Summary 


The  features  of  the  cap  model  and  the  relevant  equations  were  documented 
in  this  chapter.  In  addition,  the  steps  required  to  implement  the  cap  model 
into  the  F£  code  were  summarized.  The  implementation  of  the  cap  model 
was  verified  by  comparing  the  output  from  the  FE  code  and  a  driver  for  the 
cap  model. 
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4  Equations  of  State  for  Air, 
Water,  and  Solids 


Introduction 

The  materials  of  interest  to  this  investigation  include  partially-saturated 
soil-  and  rock-like  geomaterials  and  man-made  concretes,  grouts,  and  grout- 
cretes.  As  outlined  in  Chapter  2,  the  equations  developed  from  the  Biot 
theory  require  an  expression  for  the  bulk  modulus  of  the  pore  fluid  and  the 
grain  solids.  To  determine  the  bulk  modulus  of  the  pore  fluid,  the  concept  of 
a  homogeneous  pore  fluid  will  be  adopted  to  treat  partially-saturated  materials. 
This  investigation  will  assume  that  the  liquid  within  the  pore  space  is  water 
and  the  gas  within  the  pore  space  is  air.  Thus,  the  pore  fluid  will  be  regarded 
as  a  compressible  mixture  of  air  and  water.  Based  on  the  equations  of  state 
(EOS)  for  air  and  water,  we  will  develop  equations  for  the  bulk  modulus  of 
this  air-water  mixture.  The  grain  solids  will  be  treated  as  either  linear  or 
nonlinear  elastic  materials  or  as  nonlinear  hysteretic  materials:  each  method 
for  calculating  the  bulk  modulus  of  the  grain  solids  will  be  described. 


Equation  of  State  for  Water 

Over  the  pressure  range  of  interest  to  this  investigation,  i.e.,  0  to  600 
MPa.  water  has  a  finite  compressibility  and  should  be  treated  as  a  nonlinear 
elastic  conpressible  material.  The  compressibility  of  water  is  depicted  in 
Figure  4.1  as  a  plot  of  pressure  versus  volume  strain.  The  reader  should  note 
that  at  600  MPa  the  volume  strain  of  water  is  nominally  14  percent.  The  bulk 
modulus  or  con^ressibility  of  water  was  evaluated  from  the  Walker-Stemberg 
EOS  for  water  (Walker  ami  Sternberg  196S),  which  is  valid  for  pressure 
levels  of  up  to  SO  GPa.  The  EOS  expresses  the  water  pressure  as  an  analyti¬ 
cal  function  of  the  density  and  the  internal  energy  of  the  water.  For  this 
investigation,  the  eniergy  dependent  terms  in  the  EOS  were  not  included  due  to 
the  assumed  quasi-static  aiKl  isothermal  nature  of  the  intended  calculations. 
Without  the  energy  tenns,  the  EOS  is  expressed  as 
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Figure  4.1 .  Pressure  versus  volume  strein  response  of  water 


i 


P  ‘  Pfl  *  *  P^h  *  P^fA  -  Pq 

where  P  is  the  pressure  in  the  water,  p  is  the  density  of  the  water,  /,  are 
material  constants  and  Pq  is  the  initial  pressure.  If  we  define  volumetric 
strain  as 


then  the  bulk  modulus  of  water  may  be  expressed  as 

IL  ^  eLil  4.3 

PQ  ^P 

Substituting  Equation  4.1  into  Equation  4.3,  one  obtains  the  Fmal  expression 
for  the  bulk  modulus  of  water  as  a  fimction  of  density 

•  -  (pV,  •  3pV2  •  5pV,  *  7,V4  ) 

Po 

In  the  FE  program,  pressure  is  the  known  quantity,  not  density.  Since  the 
EOS  expresses  pressure  and  bulk  nradulus  as  a  function  of  density,  Newton's 
method  was  used  to  calculate  the  density  for  a  given  pressure,  then  the  bulk 
modulus  was  calculated. 
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Air-Water  Compressibility 

Background 

The  concept  of  a  homogeneous  pore  fluid  (HPF)  was  first  introduced  by 
Chang  and  Duncan  (1977).  In  using  their  concept,  one  assumes  that  a 
three-phase  material  containing  air,  water,  and  solids  may  be  replaced  with  a 
two-phase  material  containing  a  compressible  pore  fluid  and  solids.  A  partial¬ 
ly-saturated  material  is  transformed  into  a  fully-sanirated  material  with  a  HPF. 
Effective  stress  is  calculated  in  the  same  manner  as  for  a  fully-saturated  mate¬ 
rial,  and  the  modulus  of  the  pore  fluid  is  calculated  based  on  the  modulus  or 
compressibility  of  an  air-water  mixture.  The  concept  is  applicable  to  materi¬ 
als  with  saturation  levels  greater  than  8S  percent.  At  these  levels  of  satura¬ 
tion,  the  air  should  be  in  the  form  of  occluded  bubbles  uniformly  distributed 
throughout  the  water,  and  the  air  and  water  pressures  should  be  identical.  At 
lower  saturation  levels,  one  can  not  guarantee  that  the  air  and  water  pressures 
will  be  the  same. 

The  compressibility  of  air-water  mixtures  has  been  studied  by  several 
investigators.  Bishop  and  Eldin  (19S0)  examined  non-zero  total-stress  friction 
angles  measured  during  undrained  shear  tests.  They  attributed  the  observed 
behavior  to  incomplete  saturation  of  test  specimens.  Using  Boyle’s  and 
Hemy’s  Laws,  they  developed  expressions  for  the  compressibility  of  an  air- 
water  mixture  without  accounting  for  surface  tension  effects. 

Schuurman  (1966)  reviewed  the  work  of  previous  investigators  and  con¬ 
cluded  that  surface  tension  effects  must  be  included  in  an  air-water  compress¬ 
ibility  formulation.  Schuurman  claimed  to  be  the  first  to  attempt  such  a  for¬ 
mulation.  Schuurman  assumed  that  at  saturation  levels  greater  than 
85  percent,  the  air  existed  in  the  form  of  bubbles.  However,  to  account  for 
surface  tension,  the  radius  of  the  air  bubbles  was  required,  yet  little  if  any 
experimental  data  was  available  to  provide  this  necessary  information. 
Schuurman’ s  formulation  also  differed  from  that  of  Bishop  in  that  he  wrote  his 
expressions  in  terms  of  the  current  volume  of  air  as  opposed  to  the  original 
volume,  and  he  assumed  the  water  was  incompressible. 

Fredlund  (1976)  also  developed  an  expression  for  the  compressibility  of  an 
air-water  mixnire  using  a  formulation  in  which  the  water  had  a  finite  com¬ 
pressibility.  He  accounted  for  surface  tension  in  a  manner  that  did  not  require 
a  knowledge  of  air  bubble  sizes  by  using  a  parameter  for  air-water  pressures 
similar  to  Skempton’s  B  parameter,  which  could  be  evaluated  experimentally. 
Fredlund  also  interpreted  the  mixture  volume  in  the  expression  for  compress¬ 
ibility  as  the  volume  of  water  plus  free  air  as  conq^ared  to  water  plus  total  air. 

Chang  and  Duncan  (1977)  based  their  expressions  for  the  compressibility 
of  an  air-water  mixture  on  the  equations  of  Schuurman.  Like  Schuurman, 
they  included  surface  tension  effects  in  their  formulation  and  assumed  the 
water  was  incompressible. 
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Alonso  and  Lloret  (1982)  reviewed  the  work  of  previous  investigators, 
compared  the  conqjressibility  curves  of  each,  and  formulated  their  own  ex¬ 
pressions  for  the  compressibility  of  an  air-water  mixture.  They  assumed  a 
finite  compressibility  for  water  and  accounted  for  surface  tension  in  the  same 
manner  as  Fredlund. 

In  summary,  there  are  significant  differences  in  the  equations  developed  by 
several  investigators  for  air-water  mixtures.  For  this  reason,  equations  for  the 
con^ressibility  of  an  air-water  mixture  will  be  developed  in  this  chapter. 

Prior  to  developing  the  equations,  a  brief  description  of  the  appropriate  physi¬ 
cal  laws  will  be  provided. 


Boyle's  and  Henry's  laws 

Boyle’s  and  Henry’s  Laws  will  be  used  in  developing  equations  for  the 
compressibility  of  an  air-water  mixture.  These  laws  are  defined  and  explained 
lor  purposes  of  completeness.  Boyle’s  Law  states  that  "at  a  constant  tempera¬ 
ture,  the  volume  of  a  given  quantity  of  any  gas  varies  inversely  as  the 
pressure  to  which  the  gas  is  subjected"  (CRC  Handbook  1980). 

Air  dissolves  in  water  according  to  Henry’s  Law,  which  states  that  "the 
weight  of  gas  dissolved  in  a  fixed  quantity  of  liquid,  at  constant  temperature, 
is  directly  proportional  to  the  pressure  of  the  gas  above  the  solution"  (Fredlu¬ 
nd  1976).  Fredlund  (1976)  explains  that  the  structure  of  water  molecules  pro¬ 
duces  a  "porosity"  within  the  water  of  approximately  2  percent  by  volume. 
This  porosity  can  be  filled  by  a  gas  such  as  air,  i.e.,  air  dissolves  in  water  by 
filling  this  pore  space  (see  Table  4.1). 

Fredlund  (1976)  provides  a  simple  Table  4.1 . 
analogy  to  understand  the  congress-  Solubility  of  Air  in  Water 
ibility  of  an  air- water  mixture. 

Consider  a  test  vessel  made  of  a 
cylinder  and  piston.  At  the  base  of 
the  cylinder  is  a  porous  stone  having  a 
porosity  of  2  percent;  the  porous 
stone  simulates  the  behavior  of  the 
water.  The  piston  is  initially  posi¬ 
tioned  some  distance  above  the  stone 
with  air  filling  the  space  in  between. 

An  imaginary  valve  at  the  surface  of 
the  porous  stone  controls  the  move¬ 
ment  of  air  into  the  stone.  The  air  in 
the  porous  stone  simulates  the  air  dis¬ 
solved  in  water.  If  the  valve  is  closed 
and  the  piston  moves  down  into  the  cylinder,  the  air  above  the  stone  com¬ 
presses  following  Boyle’s  Law.  If  the  valve  is  opened,  some  of  the  air  will 
diffuse  into  the  porous  stone  following  Henry’s  Law.  This  process  will  con¬ 
tinue  until  all  of  the  air  passes  into  the  porous  stone.  When  the  piston  con- 


I  Tamperatura 

Dagraas  C 

H«nry'* 

Constant 

0 

0.0291B 

4 

0.02632 

10 

0.022B4 

15 

0.02055 

1 

0  01B6B 

1 

0.01 70B 

1  30 

0,01564 

1  from  Frtdiund  (1976)  | 

Chapiar  4  Equations  of  Staia  tor  Air,  Water,  and  Solids 


tacts  the  porous  stone,  there  is  a  discontimiity  in  the  compressibility  of  the 
system;  the  compressibility  jumps  immediately  to  that  of  water.  The  level  of 
saturation  within  an  air-water  mixture  must  be  evaluated  to  determine  the 
discontinuity  point. 


Derivation  of  equations 

The  following  assun^tions  were  made  for  this  analysis.  We  will  assume 
initial  saturation  levels  are  greater  than  85  percent,  which  implies  that  all  air 
bubbles  are  occluded.  Surface  tension  effects  will  be  neglected,  which  allows 
us  to  assert  that  the  air  bubbles  within  the  water  will  be  at  the  same  pressure 
as  the  water.  The  air  is  soluble  in  water  and  observes  Henry’s  Law,  and  the 
rate  of  increase  in  pore  water  pressure  from  any  simulation  is  slower  than  the 
rate  of  diffusion  of  air  in  water.  Finally,  prior  to  full  saturation,  the  com¬ 
pressibility  or  bulk  modulus  of  water  is  a  constant.  We  will  first  develop  the 
equations  for  an  air-water  system  with  a  rigid  porous  skeleton,  then  one  with 
a  compressible  porous  skeleton. 

The  following  terms  are  used  in  the  derivation  of  the  compressibility  of  an 
air- water  mixture.  Let 

V  denote  the  total  volume  of  air  and  water, 

V^,  the  volume  of  the  void  space, 
the  volume  of  water, 
the  volume  of  dissolved  air, 
the  total  volume  of  air, 

VJ  the  volume  of  free  air,  which  is  equal  to  -  V^, 
the  pore  water  pressure, 
the  pore  air  pressure  and 
H  Henry’s  constant. 

A  subscripted  "o"  is  used  to  indicate  an  initial  value.  The  total  mass  of  air 
and  water  remains  constant.  Substituting  expressions  for  porosity  (n)  and  sat¬ 
uration  (S),  the  initial  volumes  of  water  and  free  air  may  be 
expressed  as 

‘•5 


-S„)V.„  =  n„V„(l  -S„) 


Using  Henry’s  Law  and  Equation  4.5,  the  initial  volume  of  dissolved  air  may 
be  expressed  as 

Vj,  ■  K,,  4,7 


The  sum  of  Equations  4.6  and  4.7  yields  an  expression  for  the  initial  toul 
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volume  of  air  in  the  system 

4.8 

Expressions  for  the  compressibility  of  water  and  an  air-water  mixture  may  be 
written  as 


C 


y.  dp^ 


C  I  dV,  ^  dV^ 

v  ' ,  V  dP,  dP, 
a  w  I 


respectively.  Substituting  Equation  4.9  into  Equation  4.10  one  obtains 


=  - 
ftl 


1  dV^ 

^  -  V  r 

V  ♦  V 

a  w 


We  will  now  use  Boyle’s  Law  to  develop  an  expression  for  the  derivative 
in  Equation  4.11.  Boyle’s  Law  may  be  written  as 


V  P  =  V  P 

'  a  ao  '  ao 


If  we  assume  the  pore  and  air  pressure  are  equal  and  we  can  write 

the  following 


from  which  we  write 


Pao  P.o  Vg  *  Vd 

Pg  Pw  V  *  Vj 

go  d 


y.n  *  y^ 

go  a 


V'.  ^  yd 


V  P 

ao  *  ao 


Substituting  the  laner  expression  in  Equation  4  14  into  Equation  4. 1 1 ,  one 
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obtains 


Va  ^  K  i 


which  is  an  expression  for  the  compressibility  of  an  air- water  mixture.  By 
judiciously  substituting  for  the  volume  terms  in  Equation  4.15,  we  will  devel¬ 
op  a  final  expression  for  the  compressibility  of  the  mixture. 

By  combining  Boyle’s  Law  (  Equation  4.12)  and  Equation  4.8,  we  may 
write  an  expression  for  the  current  total  volume  of  air 

K.  =  4.16 

“a 

The  current  volume  of  water  may  be  expressed  as 

V  =  V  ( 1  +  C  6P)  4.17 

and,  after  substituting  for  as 

The  current  volume  of  dissolved  air,  which  is  a  ftmction  of  Henry’s  Law  and 
the  current  volume  of  water,  is  written  as 

n,V„S,H  il  .  C^5P)  4.19 

Subtracting  Equation  4.19  ft’om  Equation  4.16,  one  obtains  an  expression  for 
the  current  volume  of  free  air 


I'a  -  n^VA^il  -S„*S„H)  -  C^6P) 


Adding  Equations  4.20  and  4.18,  one  obtains 


Va  *  I'k  =  "oVo  *  S„(  1 -//)(!  .  C^5P) 

4.21 


which  will  eventually  be  substituted  back  into  Equation  4.15.  Combining 
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Equations  4.8  and  4.16,  one  may  write 


and,  by  multiplying  Equation  4.18  by  the  compressibility  of  water  and 
dropping  the  higher  order  terms,  one  obtains 


Substituting  Equations  4.21,  4.22,  and  4.23  into  Equation  4.  IS  yields  the  final 
expression  for  the  compressibility  of  an  air-water  mixture 


C„-  \^{i-  S,H)  *  S^(l-  H){1-  C^bP) 

‘a 


In  a  similar  manner,  an  expression  for  the  level  of  sahiration  may  be  devel¬ 
oped  and  written  as 


V  +  V 


S,(1*C,6P) 


{\-S^*S^H)  *  S^{\-  H){\*  C^bP) 


When  the  porous  skeleton  is  compressible,  the  current  void  volume  may  be 
expressed  as 


«**> 


where  is  the  initial  toul  volume  of  voids  and  solids  and  t  **  is  the  effec¬ 
tive  volumetric  strain.  Substituting  the  above  and  Equation  4.18  into  the  first 
expression  in  Equation  4.25,  one  obuins 
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n^S,(l-C^dP) 
%  ~  ^kk 


4.27 


S  = 


W 


V  *  V 

w 


which  is  an  expression  for  the  level  of  saturation  in  a  deforming  porous  skel¬ 
eton.  An  equation  for  the  compressibility  of  an  air-water  mixture  within  a  de¬ 
forming  porous  skeleton  may  be  obtained  by  combining  Equations  4.15,  4.22, 
4.23,  and  4.26  to  yield 


Cm  = 


%  ~ 


p 

(1 


Sofi) 


4.28 


In  the  process  of  a  calculation,  one  must  first  evaluate  Equation  4.27.  If  • 

the  level  of  saturation  is  less  than  one.  Equation  4.28  is  used  to  calculated  the 
bulk  modulus  of  the  pore  fluid.  If  the  level  of  samration  is  equal  to  one,  the 
bulk  modulus  of  the  pore  fluid  is  calculated  from  the  EOS  of  water,  i.e.. 

Equation  4.4. 

To  illustrate  the  response  of  a  partially-samrated  material  to  an  ^plied  * 

loading,  an  example  calculation  was  conducted  and  the  output  graphically 


Figure  4.2.  Pressure-volume  response  of  partially-saturated  material 


presented  in  Figure  4.2.  The  simulated  material  has  a  Young’s  modulus  of 

1800  MPa.  a  bulk  modulus  of  1000  MPa.  a  total  porosity  of  20  percent,  a  • 
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A 


saturation  level  of  90  percent,  and  an  air  porosity  of  2  percent.  The  material 
was  loaded  under  undrained  uniaxial  strain  boundary  conditions.  At  volume 
strains  less  than  approximately  2  percent,  the  generated  pore  pressures  are 
negligible,  and  the  loading  bulk  modulus  is  equal  to  the  skeletal  bulk  modulus 
of  1000  MPa.  At  these  strains  levels,  the  material  loads  as  if  it  were  fully 
drained.  At  a  volume  strain  of  2  percent,  all  of  the  air  porosity  is  eliminated, 
and  the  pore  fluid  becomes  fully  saturated  and  begins  to  carry  a  major  portion 
of  the  applied  stress.  At  these  strain  levels  and  above,  the  material  loads  as  a 
fully-saturated  material.  In  addition,  the  pressure-volume  response  is  nonlin¬ 
ear  due  to  the  nonlinear  namre  of  the  water. 


Equation  of  State  for  Solids 

Three  methods  for  calculating  the  bulk  modulus  of  the  grain  solids  were 
inplemented  into  the  FE  program.  The  first  method  assumes  the  grain  solids 
are  linear  elastic  materials.  The  second  method  uses  an  analytical  EOS  and 
treats  the  solids  as  a  nonlinear  elastic  material.  The  third  method  uses  a 
simple  model  to  simulate  the  nonlinear  hysteretic  material  behavior  of  the 
graiiis. 

The  first  method  is  self  explanatory;  the  program  sinply  uses  a  constant 
bulk  modulus  value  for  the  entire  calculation.  In  the  second  method,  an  ana¬ 
lytic  relationship  between  pressure  and  compression  is  developed  for  each 
material.  Compression  is  defined  as 

u  =  .  *  4.29 

I  -  t 

where  t  is  the  Cauchy  or  engineering  strain.  Using  solid  carbonate  as  an 
example  material,  the  pressure-compression  relationship  is  linear  below 
1 .2  GPa  and  may  be  written  as 

Pg  "  0.7  n  4.30 

where  P,  is  the  grain  pressure.  The  bulk  modulus  for  carbonate  may  then  be 

A 

written  as 

K.  =  (  1  -  M  )*  =07(1-^)*  4.31 

•  dn 

A  plot  of  pressure  versus  volumetric  strain  for  carbonate  is  plotted  in 
Figure  4.3  Other  materials  may  be  simulated  in  an  analogous  manner. 

The  third  model,  which  simulates  nonlinear  hysteretic  material  behavior, 
uses  tabulated  curves  that  describe  the  loading  and  unloading  pressure-volume 
response  of  the  grains.  This  model  is  based  on  the  work  of  Meier  ( 1986), 
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Figure  4.3.  Pressure  versus  volume  strein  respor^se  of  cerbor^ate 


who  used  a  similar  model  fur  one-dimensional,  plane-wave  ground  shock 
calculations. 

During  virgin  loading,  the  volume  strain  is  computed  from  the  current 
value  of  pressure  through  linear  interpolation  of  the  tabulated  loading  curve. 

A  similar  process  is  employed  when  unloading  occurs  from  pressure  levels  at 
or  above  the  lockup  point  using  the  tabulated  unloading  curve.  When  unload¬ 
ing  takes  place  from  pressure  levels  below  the  lockup  point,  a  scaling  process 
must  be  applied  to  the  tabulated  unloading  curve.  In  this  scaling  process,  let 
and  („  represent  the  peak  pressure  and  peak  volume  strain,  respectively, 
from  which  unloading  commences  (see  Figure  4.4).  The  pressure  and  voIuhk 
strain  scaling  factors  are  calculated  as 


and 

4.33 

where  /*/  is  the  pressure  at  lock  up.  F,  is  the  tension  cutoff  pressure,  and  o  is 
an  empirically  determined  calibration  constant  with  values  ranging  between 
:mro  and  unity.  Knowing  the  unloading  pressure  (/*,),  the  recovered  pressure 

(A  P).  which  is  the  difference  between  the  pressure  »:  lock  up  and  the  value 
of  pressure  on  the  tabulated  unloading  curve,  is  calculated  in  the  following 
manner 


60 


Chapt*«  4  EquatioTM  ot  Siai*  lot  Ait.  Watot.  and  SoMt 


Figure  4.4.  Nonlinaer-hytteretic  model 


4.34 


The  recovered  strain  ( A « )  is  conned  through  linear  interpolation  of  the 
tabulated  unloading  curve.  The  unloading  volume  strain  (c ,)  is  calculated  by 
subtracting  the  scaled  value  of  recovered  strain  horn  the  peak  strain 


The  bulk  modulus  on  this  unloading  curve  is  calculated  as 


K.ifjfp) 


K, 


1  -  a 


AM 


where  >  AP I  At.  Reloading  occun  along  a  line  passing  between  the 
last  unloading  pressure-volume  strain  point  and  the  point  P„,f„. 
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Summary 


The  algorithms  required  to  numerically  simulate  the  behavior  of  the  three 
primary  constituents  of  geomaterials,  air,  water  and  solids,  were  documented 
in  this  chapter.  When  these  algorithms  are  combined  with  an  appropriate 
skeletal  model  into  the  F£  formulation  of  Biot’s  theory,  the  multikilobar 
response  of  any  geomaterial  may  be  calculated. 
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5  Features  and  Verification  of 
FE  Program 


Introduction 

The  objectives  of  this  ct^pter  are  to  (i)  describe  features  in  the  FE  pro¬ 
gram  that  have  not  already  been  introduced  and  (2)  present  solutions  from 
several  verification  problems  as  proof  that  the  program  works  correctly. 
Several  features  of  ^e  FE  program  have  already  been  introduced.  In  Ch^ 
ter  2,  the  benefits  gained  ftom  effective  stress  simulations  of  multi-kilobar 
material  behavior  and  the  available  material  models  were  described.  The  cap 
model  was  documented  in  Qiapier  3  and  the  equations  of  state  of  air.  water, 
and  solids  were  described  in  Chapter  4. 

The  following  features  will  be  presented  in  this  chapter.  A  restart  feature 
was  implemented  into  the  FE  program  to  permit  the  simulation  of  certain 
laboratory  tests.  For  exan^le,  a  consolidated  tmdrained  triaxial  compression 
test  wherein  the  consolidation  phase  has  drained  boundary  conditions  and  the 
shear  phase  has  undrained  boundary  conditions  requires  changing  fluid  flux 
boundary  conditions.  A  brief  summary  of  the  postprocessing  procedures  will 
also  be  presented.  These  features  are  described  in  the  next  section. 

The  final  sections  of  this  chapter  document  solutions  from  several 
verification  problems.  For  each  problem,  the  FE  solution  is  compared  to 
a'  aitable  cl(^  form  or  analytic  solutions.  These  verification  problems 
e  iblish  the  FE  program's  ability  to  correctly  solve  a  variety  of  initial  and 
boundary  value  problems. 


Additional  Features  of  FE  Program 

Rastart  fMtura 

The  rt^tait  feature  was  itr^lonented  for  the  purpose  of  allowing  the  user  to 
change  the  boundary  conditions  at  a  preselected  time  in  the  calculation.  A 
K(,/BX/STX  test  (Kronyms  defined  subsequently)  is  an  example  of  a  labora¬ 
tory  test  with  changing  boundary  conditions.  This  test  is  conducted  by 
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loading  a  cylindrical  specimen  to  a  desired  mean  normal  stress  level  under  K,, 
or  uniaxial  strain  boundary  conditions,  unloading  to  a  desired  mean  normal 
stress  level  under  constant  axial  strain  (BX)  boundary  conditions,  and  then 
conducting  a  constant  radial  stress  triaxial  compression  (SIX)  test  at  yet 
another  mean  normal  stress  level.  The  loading  and  the  BX  unloading 
phases  may  be  numerically  simulated  with  displacement  controlled  boundary 
conditions.  However,  to  realistically  attempt  to  simulate  the  STX  phase,  the 
user  should  apply  stress  controlled  boundary  conditions.  The  restart  feature 
implemented  into  JAM  allows  the  user  to  perform  this  calculation  in  a  simple 
manner. 


Postprocessing 

In  many  instances,  one  would  like  to  plot  FE  results  at  a  single  location, 
e.g.,  at  the  nodal  points.  Many  postprocessing  FE  software  packages  require 
stress  and  strain  values  at  the  nodes  rather  than  at  the  Gauss  integration 
points.  A  procedure  was  implemented  in  the  FE  program  JAM  to  extrapolate 
and  smooth  Gauss  point  data  to  the  element  vertices,  i.e.,  comer  nodes. 
Values  at  the  midside  nodes  were  then  calculated  from  the  values  at  the 
appropriate  comer  nodes. 

The  implemented  smoothing  procedure  was  developed  and  described  by 
Hinton.  Scott,  and  Ricketts  (1975)  and  Hinton  and  Campbell  (1974).  The 
procedure  is  simple  and  straightforward.  The  smoothed  stresses  at  the  nodes 
may  be  calculated  ft'om  the  expression 
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where  the  a ,  are  the  smoothed  stresses,  o j.  a jj.o m  and  0 are  the  stresses 

at  the  integration  points,  and  a  ^  \  ^  .  b  •  ~  1  and  c  =  1  - 

At  a  given  comer  node,  smoothed  values  from  adjacent  elements  are  averaged 
to  yield  a  single  value. 


Plane  and  Axisymmetric  Verification  Problems 

Five  problems  with  plane  or  axisymmetric  geometries  were  solved  with 
JAM  to  test  and  verify  that  the  material  models,  the  eight-node  quadrilateral 
element,  and  numerous  other  algorithms  were  correctly  implemented  in  the  FE 
program.  For  each  of  tlK  five  problems,  selected  output  from  JAM  are 
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Verification  Problem  1  exercises  a  single  element  under  distributed  normal 
and  shear  loads  of  1000/length  as  shown  in  Figure  5.1.  The  element  sim¬ 
ulates  an  isotropic  linear  elastic  material  having  a  Young’s  modulus  of  30  x 

10  *  and  a  Poisson’s  ratio  0.3.  The  following  boundary  conditions  were  • 

imposed: 

«  0  at  point  A  and  <  0  at  points  B  and  C 


Tabl«5.1. 

Result*  from  Vsrification  Probism  1 


Table  5. 1  compares  stress  and  strain  states  calculated  from  the  FE  and 
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analytic  solutions  (Hibbitt,  Karlsson  and  Sorensen  1989)  for  this  problem 
under  plane  strain  and  plane  stress  boundary  conditions.  This  problem 
exercises  the  elastic  constitutive  model,  verifies  that  the  eight-node  quadratic 
element  accurately  models  constant  strain  states,  and  also  checks  that 
distributed  loads  are  correctly  simulated.  The  results  from  JAM  match  the 
analytic  solution  exactly.  Verification  Problem  1  was  also  successfully  solved 
using  the  Cap  model  algorithm. 


Verification  Problem  2  exercises  a  single  axisymmetric  element  under 
distributed  normal  loads  of  1000/arca  as  shown  in  Figure  5.2,  The  element 
simulates  an  isotropic  linear  elastic  material  having  a  Young’s  modulus  of 
30  X  10  ^  and  a  Poisson's  ratio  0.3.  The  following  boundary  conditions  were 
imposed: 

u.  =  0  at  points  A,  B,  and  C 


Table  5.2  compares  stress  and  strain  sutes  calculated  from  the  FE  and 
analytic  solutions  (Hibbitt,  Karlsson  and  Sorensen  1989)  under  the  imposed 
axisymmetric  boundary  conditions.  Like  Verification  Problem  1,  problem  2 
exercises  the  elastic  constitutive  model,  verifies  that  the  eight-node  quadratic 
element  accurately  models  constant  strain  sutes,  and  also  checks  that 
distributed  loads  are  correctly  simulated.  The  results  from  JAM  match  the 
analytic  solution  exactly. 

Plane  and  axisymmetric  patch  tests  were  employed  in  Verification 
Problems  3  and  4.  In  the  patch  test,  nodal  point  displacements  are  an^lied  to 
a  patch  of  elements  such  that  a  constant  swe  of  strain  exists  throughout  the 
mah.  In  Verification  Prt^lem  3,  the  elemems  simulate  an  isotropic  linear 
elastic  material  having  a  Young's  modulus  of  30  x  lO^  and  a  Poisson's  ratio 
0.3.  The  impened  displacement  boundary  conditions  were  appliol  to  the  patch 
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Table  5.3. 

Results  from  Verification  Problem  3 


Strata 

or 

Strain 

Plana  Strain 

Plano  Strata 

JAM 

Analytic 

JAM 

Analytic 

1600. 

1600. 

1333. 

1333. 

1600. 

1600. 

1333. 

1333. 

400. 

400. 

400. 

400. 

800. 

800. 

0. 

0. 

1.x  10-3 

1.x  10-3 

1 .  X  1 0-3 

1 .  X  1 0'3 

'y 

1.  X  10-3 

1.x  10'^ 

1.x  10'^ 

1.x10‘3 

1.  x10'3 

1.x  10-3 

1 .  X  1 0-3 

1.x  10-3 

simulatine  an  isotropic  linear  elastic  material  having  a  Young’s  modulus  of 
30  X  10”  and  a  Poisson’s  ratio  0.3.  The  imposed  displacement  boundary 
conditions  were  applied  to  the  patch  of  elements  shown  in  Figure  5.4  and 
were  calculated  as; 


m  m 

(r  -  1000)  >  1 


xlO'^  and  «.  = 
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(r-  1000)' 
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-3 


Table  5.4  con^sares  stress  and  strain  sutes  calculated  from  the  FE  and 
analytic  solutions  (Hibbitt.  Karlsson  and  Sorensen  1989)  for  this  problem. 
The  results  from  JAM  match  the  analytic  solution  exactly.  Verification 
Problem  4  was  also  successfully  solved  using  the  Cap  model  algorithm. 
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Table  5.4. 

Results  from  Verification  Problem  4 


Strasa 

or 

Strain 

Axityinmatric 

JAM 

An^ytic 

5.769x10^ 

b.769x  10^ 

<^1 

5.769x10^ 

5.769X  10^ 

1.154X  10^ 

1 , 1 54  X  1 0^ 

3.462  X  10^ 

3.462  X  10^ 

<f 

1.x  10'^ 

1 .  X  1 0  3 

1 .  X  1 0-3 

1.x  10'^ 

1.x  10-3 

l.xlO'3 

0. 

0. 

Verification  Problem  S  is  a  plane  strain  simulation  of  a  thick  wall  cylinder 
subjected  to  an  increasing  internal  pressure.  Due  to  the  symmetry  of  the 
problem,  a  quarter  grid  was  used  in  the  calculation;  the  problem  geometry  and 
FE  mesh  are  shown  in  Figure  S.S.  The  material  was  modeled  with  the 
following  properties,  a  Young’s  modulus  of  21000,  a  Poisson’s  ratio  of  0.3,  a 
yield  stress  of  24  and  a  linear  hardening  modulus  of  0.  When  a  cylinder  with 
these  properties  is  subjected  to  an  increasing  internal  pressure  above  10.4,  an 
elastic-plastic  boundary  moves  through  the  cylinder;  on  the  external  side  of 
the  boundary,  all  of  the  strains  are  elastic,  and  on  the  interior  side,  the  strains 
are  elastic-plastic.  Table  S.S  compares  stresses  calculated  from  the  FE  and 
analytic  solutions  (Hodge  and  White  i9S0;  Prager  and  Hodge  i9Si)  at  several 
radii  within  the  elastic  region  for  an  applied  imemal  pressure  of  18,  which 
places  the  elastic-plastic  boundary  at  a  radius  of  160.  The  computed  results 
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Table  5.5. 

Results  from  Verification  Problem  5 


Straii 

Radius 

Plana  Strain 

JAM 

Analytic 

Max. 

Principal 

Stress 

or 

163.64 

22.12 

22.09 

176.32 

20.32 

20.25 

193.64 

18.37 

18.31 

>  160 

5.33 

5.31 

>  160 

23.11 

23.03 

are  well  within  the  imposed  convergence  tolerance  of  1  percent.  The  results 
from  JAM  also  agree  with  the  results  calculated  by  Owen  and  Hinton  (1980) 
for  this  problem.  This  validates  the  plasticity  formulation  in  JAM. 

Consolidation  Problems 

To  verify  that  the  FE  progiam  could  solve  consolidation  problems,  output 
from  JAM  were  compared  to  the  results  calculated  from  closed  form  solu- 


Figure  5.6.  Pore  pressure  versus  depth  at  two  time  increments 

tions.  Boundary  conditions  and  material  properties  were  altered  to  fully  exer¬ 
cise  different  features  within  the  FE  program.  In  Figure  S.6,  depth  versus 
calculated  pore  fluid  pmsures  are  plotted  in  a  normalized  format  at  two  dif- 
ferem  time  incremems  for  a  one-dimensional  consolidation  problem  in  whidi 
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the  soil  column  was  idealized  as  an  elastic  porous  skeleton  with  an  incompres¬ 
sible  pore  fluid.  Good  agreement  is  shown  between  the  FE  results  and  the 
closed  form  solution.  A  similar  one-dimensional  problem  was  solved  with 
two  materials  having  compressible  pore  fluids,  where  the  ratios  of  pore  fluid 
modulus  to  skeletal  modulus  (N)  were  2000  (nearly  incompressible  pore  fluid) 
and  5  (highly  compressible  pore  fluid).  Results  from  the  FE  program  and  the 
closed  form  solution  (Chang  and  Duncan  1983)  are  plotted  in  Figure  5.7  as 
normalized  displacements  versus  time  factor,  i.e.,  normalized  time.  Again, 


Figure  5.7.  Displacement  versus  time  for  one-dimensional  consolidation  of 
an  idealized  elastic  material 

the  results  show  reasonable  agreement  between  the  FE  results  and  the  closed 
form  solution. 

A  two-dimensional  axisymmetric  consolidation  problem  consisting  of  a 
circular  foundation  on  a  finite  soil  layer  (Figure  S.8)  was  also  calculated.  The 
mesh  is  A  units  high  by  lOA  units  wide,  and  a  uniform  vertical  load  of  radius 
A  was  applied  to  the  top  surfaces  of  three  elements  to  simulate  the  foundation 
loads.  The  following  boundary  conditions  were  invoked  for  this  problem. 

The  vertical  edges  of  the  mesh  (A-D  and  B-C)  were  constrained  in  the  radial 
direction,  the  bottom  edge  of  the  mesh  (C-D)  was  constrained  in  the  vertical 
direction,  the  top  surface  (A-B)  was  free  draining,  and  no  flow  conditions 
were  applied  to  the  three  remaining  surfaces  (B-C.  C-D.  and  A-D).  The  cal¬ 
culated  senlements  (in  dimensionless  format)  from  JAM  (solid  circles)  are 
compared  in  Figure  5.9  to  settlements  calculated  from  the  analytical  solution 
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Figure  5.8.  Mesh  geometry  for  Axisymmetric  Consolidation  Problem 


Figure  5.9  Displacement  vs  Time  from  20  Consolidation  Problem 

(solid  line).  Again,  there  is  excelleffi  agreement  between  the  calculuion  and 
the  analytical  solution. 
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Cryer  Problem 

A  numerical  simulation  of  Cryer ’s  problem  (Cryer  1%3)  was  conducted  as 
an  additional  verification  test  of  the  FE  program.  Cryer  developed  a  closed 


Figure  5.10.  Mesh  geometry  for  Cryer  Problem 


Figure  5.11 .  Dimensionless  pore  pressure  response  for  Cryer's  problem 


form  solution  to  the  problem  of  i  sphere  of  elastic  purous  matcnal  loaded  on 
the  surface  by  a  constaiu  uniform  pressure  and  having  drained  boundary  con¬ 
ditions.  For  vaJiKS  of  Poisson's  ra  lO  less  than  0.5.  pore  pressure  at  the 
center  of  the  sphere  increases  to  stress  levels  greater  than  the  externally 
ai^lied  pressure  and  then  dissipates.  The  greatest  increase  in  pore  pressure 
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Utilizing  the  symmetry  of  the  problem,  the  soil  sphere  was  represented  by 
the  mesh  depicted  in  Figure  5.10  and  calculated  as  an  elastic  axisymmetric 
problem.  A  unit  pressure  was  placed  on  the  boundary  during  the  first  incre¬ 
ment  of  loading  and  held  constant  thereafter.  Figure  5.11  and  Figure  5.12 
compare  the  results  from  JAM  with  the  closed  form  solution  for  a  value  of 
Poisson’s  ratio  equal  to  0.  Figure  5.11  is  a  plot  of  dimensionless  pore 
pressure  at  the  center  of  the  sphere  versus  the  square  root  of  dimensionless 
time;  Figure  5.12  is  a  plot  of  dimensionless  displacement  of  th“  outer  surface 
of  the  sphere  versus  the  square  root  of  dimensionless  time.  The  comparison 
between  the  FE  calculation  and  the  closed  form  solution  is  very  good.  Pore 
pressure  contours  at  a  dimensionless  time  of  approximately  0.04  are  plotted  in 
Figure  5.13.  At  this  early  time,  a  significant  portion  of  the  sphere  has  pore 
pressures  greater  than  unity. 


Summary 

In  this  chapter,  the  restart  and  post-processing  features  of  the  FE  program 
were  describwl.  The  documented  verification  problems  indicate  that  the  FE 
program  correctly  calculates  one-  and  two-dimensional  consolidation  problems 
and  elastic  and  elastic-plastic  boundary  value  problems.  Although  the  success¬ 
ful  calculation  of  these  verification  problems  docs  not  certify  the  FE  program 
is  error  free,  they  should  increase  the  confidence  of  the  end  user. 
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6  Numerical  Simulations 


Introduction 

Numerical  simulations  of  limestone  behavior  under  drained  and  undrained 
boundary  conditions  are  presented  in  this  chapter.  The  ability  of  the 
14-parameter  c^  model  to  simulate  the  basic  drained  behavior  of  limestone  is 
demonstrated  by  conq)aring  calculated  responses  of  hydrostatic,  uniaxial 
strain,  and  triaxial  con^ression  loadings  with  measured  or  recommended  lime¬ 
stone  responses.  In  a  similar  manner,  the  ability  of  the  FE  code  to  calculate 
the  undrained  behavior  of  limestone  is  demonstrated  by  comparing  calculated 
responses  of  uniaxial  strain  loadings  with  reconunended  limestone  responses. 
Finally,  some  example  calculations  are  documented  that  demonstrate  the  utility 
of  the  FE  code  in  analyzing  laboratory  test  specimen  conditions. 


Salem  Limestone 

The  limestoiM  simulated  in  this  chapter  is  commonly  referred  to  as  Salem, 
Bedford,  or  Indiana  limestone.  It  was  extraaed  from  the  Salem  formation 
near  the  community  of  Bedford,  Indiana.  Mechanical  property  tests  were 
conducted  on  intact  specimens  of  i3.S-percent  porosity  Salem  limestone  by  the 
Vermont  office  of  Applied  Research  Associates.  These  mechanical  property 
tests  included  drained  and  undrained  (with  pore  pressure  measurements) 
hydrostatic  loading  tests.  or  uniaxial  strain  tests,  triaxial  compression  tests, 
and  strain  path  tests.  Laboratory  test  diua  and  recmnmended  mateiial  proper¬ 
ties  were  obuined  from  Blouin  and  Chitty  (1988a.  1988b)  and  Zelasko  (1991). 


Simulations 

Procass 

Prior  to  numerically  simulating  limestone  behavior  under  drained  or 
undrained  bmmdary  conditions,  the  skeletal  or  drained  behavior  of  Salem 
limestone  was  required.  The  14-paranieter  cap  model,  which  was  documented 
in  Chapter  3.  was  fit  to  recommoided  drained  Salem  limestone  mechanical 
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properties.  With  this  model  and  fit  inyilemented  into  JAM,  drained  single¬ 
element  boundary  value  problems  Vk^ere  conducted  to  insure  the  FE  code 
would  reproduce  the  cap  model  calculations.  Undrained  single-element 
boundary  value  problems  were  then  conducted. 


Drained  limestone  simulations 

The  following  recommended  drained  Salem  limestone  mechanical  proper¬ 
ties  were  available  for  fitting:  a  failure  envelope,  hydrostatic  load  and  unload 
behavior,  stress-strain  behavior.  Kg  pressure-volume  behavior  and  Kg 
stress  paths,  stress-strain  curves  from  triaxial  compression  tests  conduct^  at 
several  confining  stress  levels,  and  strain  path  data  along  three  different  paths. 

Typically,  a  high  fidelity  fit  of  both  the  hydrostatic  loading  and  Kg,  or 
uniaxial  strain,  responses  is  inqxissible  to  capture  with  a  relatively  simple  cap 


I _ 

Figure  6.1.  Drained  Kg  tlress-sirain  comparison 

model.  For  this  reason,  greater  emphasis  was  placed  on  fining  the  uniaxial- 
strain  stress  path  and  stress-strain  responses  and  less  emphasis  on  the  hydro¬ 
static  loading  response.  In  Figures  6  1-6.3,  the  calculated  drained  Kg  stress 
and  strain  responses  from  the  14-parameter  cap  model  arc  compared  to  the 
recommended  drained  K,  behavior.  Figure  6. 1  compares  the  drained  Kg 
stress-strain  behavior.  Figure  6.2  the  Kg  stress  paths,  and  Figure  6.3  the  Kg 
pressure-volume  responses.  The  quality  of  the  fits  are  very  good.  To  make 
these  fits,  one  must  compromise  between  fining  the  Kg  stress  path  and  the  Kg 
stress-strain  behavior.  The  model  Kg  stress-strain  response  breaks  over  at  a 
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In  Figure  6.4,  the  calculated  drained  hydrostatic  pressure-volume  response 
of  Salem  limestone  is  con^ared  to  the  recommended  behavior.  The  quality  of 


Figure  6.4.  Drained  hydrostatic  load-unload  comparison 


tlie  fit  is  not  very  good  because  greater  emphasis  was  placed  on  fitting  the  K,, 
behavior.  Only  a  very  conqjlicated  cap  model,  with  several  tens  of  fining 
parameters,  would  fit  both  the  hydrosutic  and  behavior  of  this  material. 

Drained  triaxial  con^ression  (TXC)  tests  at  confuting  pressures  of  100  and 
400  MPa  were  also  simulated  with  the  c^  model.  The  calculated  responses 
are  ploned  as  principal  stress  difference  versus  axial  strain  and  compared  to 
actual  test  results  in  Figures  6.S  and  6.6.  The  quality  of  the  fits  is  quite  good 
considering  the  error  introduced  into  the  calculations  by  the  lack  of  fidelity  in 
the  calculated  hydrosutic  pressure- voli  me  response. 


Undrained  Hmostona  aimulations 

The  following  single-element  undrained  simulations  were  performed  using 
tt»  Walker-Stemberg  EOS  for  water  and  a  carbonate  EOS  for  the  grain  solids 
The  cap  model  fit  to  the  recommended  drained  limestone  properties  modeled 
the  skeleul  behavior  of  the  limestone. 

An  undrained  test  conducted  on  a  fully-saturated  specimen  of  Salem 
limestone  was  simulated  with  the  FE  code.  The  output  is  compared  to  the 
available  recommended  properties  as  another  method  of  verifying  the  FE 
code.  The  calculated  and  recommended  stress-strain  responses  are  compared 
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Figure  6.5.  Drained  TXC  stress-strain  comparison  at  100  MPa  confining 

pressure 


Figure  6.6.  Drained  TXC  strass-strwn  comparison  at  400  MPa  confining 

pressure 

in  a  plot  of  total  vertical  stress  versus  total  vertical  strain  (Figure  6.7).  The 
calculated  undrained  stress>strain  re^x>nse  replicates  the  recoromended 
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Figure  6.8.  Undreined  stress  path  comparison 


stress-strain  roiponse  perfmly  during  the  loading  phase,  while  the  unloading 
slightly  stiffer.  Tlw  caloUated  and  reconunended  stress  {»ths  are  companKi 
a  plot  of  princifud  stress  difference  versus  total  mean  normal  stress 
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(Figure  6.8).  The  correlation  between  the  calculated  and  the  recommended 
stress  paths  is  excellent. 

The  total  (solid),  effective  (short  dash)  and  pore  fluid  stresses  (long  dash) 
for  the  simulated  undrained  test  are  presented  in  Figure  6.9  in  the  format 


Figure  6.9.  Stresses  in  a  simulated  undrained  test 

of  stress  versus  total  volume  strain.  This  figure  illustrates  that,  even  in  a 
competent  rock  such  as  limestone,  a  significant  portion  of  the  total  applied 
stress  IS  carried  by  the  pore  fluid,  and  the  peak  effective  stress  is  only  20%  of 
the  peak  total  applied  stress. 

In  "convemional"  soil  mechanics,  water  and  the  grain  solids  are  often 
assumed  to  be  incompressible.  These  assumptions  have  significant  im^'lica- 
tions  with  regard  to  the  response  of  materials  during  undrained  loading. 

Under  undrained  or  zero  volume  change  boundary  conditions,  the  undrained 
strength  at  failure  and  the  undrained  effective  stress  path  are  unique  for  a 
given  material  with  prescribed  initial  conditions  (Lambe  and  Whitman  1%9). 
This  means  that  the  effeoive  stress  path  is  independem  of  the  applied  total 
stress  path.  A  path  of  zero  volume  change  in  ui  elastic-plastic  material 
implies  that  the  elastic  and  plastic  volume  strains  are  of  equal  magnimde  and 
opposite  sign.  To  demonstrate  that  a  unique  effective  stress  path  is  developed, 
an  undrained  triaxial  compression  (TXC)  test,  i.e. .  constant  radial  stress 
durin:  shear,  and  an  undrained  constant  mean  normal  stress  (CP)  test,  i.e., 
constant  mean  normal  stress  during  shear,  were  numerically  siimilated.  The 
following  cxmditions  existed  prior  to  the  undrained  loading  in  both  simula¬ 
tions:  (1)  the  compressibilities  of  the  water  and  the  grain  solids  were  zero; 
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(2)  the  nominal  effective  confining  stress  in  the  specimens  was  200  MPa, 
which  was  generated  by  a  drained  hydrostatic  loading;  and  (3)  the  initial  pore 
fluid  pressure  was  zero.  The  stresses  during  the  shear  loading  were  applied 
incrementally  until  the  calculation  would  not  converge  under  a  convergence 
tolerance  of  0.5  percent.  The  calculated  total  and  effective  stress  paths  for  the 
TXC  test  (solid  line)  and  the  CP  test  (dashed  line)  are  presented  in  Figure 
6. 10.  The  calculated  effective  stress  paths  from  the  TXC  and  CP  tests  are 
identical. 

Additional  undrained  calculations  were  performed  to  prove  that  the  un¬ 
drained  effective  stress  paths  are  not  unique  when  the  water  and  grain  solids 
are  compressible.  Three  undrained  TXC  tests  with  the  following  initial  condi¬ 
tions  were  simulated:  the  nominal  effective  confining  stress  was  200  MPa  and 
the  applied  back  pressures  (initial  pore  fluid  pressures)  were  0.  100,  and 
300  MPa.  The  calculated  effective  stress  paths  are  presented  in  Figure  6.11. 
These  calculations  indicate  that  as  the  applied  back  pressure  increased  from  0 
to  300  MPa,  the  corresponding  effective  stress  paths  moved  to  lower  values  of 
effective  mean  normal  stress.  This  response  can  be  explained  with  the  follow¬ 
ing  logic.  As  the  water  becomes  stiffer  with  increasing  levels  of  back  pres¬ 
sure,  «]ual  strain  increments  within  the  specimen  generate  larger  incremems 
of  pore  fluid  pressure.  Thus,  the  effective  stress  paths  move  to  the  left  in 
Figure  6.1 1. 

The  previous  sections  show  that  the  FE  code  can  accurately  simulate  both 
drained  aiKl  undrained  responses  of  .Salem  limestone  under  ideal  laboratory 
test  boundary  condition.  In  the  following  sections,  non-ideal  boundary 
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conditions  that  actually  exist  in  the  laboratory  tests  will  be  simulated. 


Test  Specimen  Simulations 

FEgrid 

Cylindrical  test  specimens  were  simulated  with  the  axisymmetric  FE  grid 
depicted  in  Figure  6.12.  A  quarter  grid,  consisting  of  144  elements  and  483 
nodes,  was  used  in  the  simulation  due  to  the  symmetric  nature  of  the  problem. 
The  specimen  end  caps  were  iiKluded  in  the  simulation  to  investigate  the 
effecu  of  end  cap  restraint  upon  both  the  stress  and  strain  conditions  within 
the  test  specimen.  A  worst  case  situation  was  simulated,  i.e.,  one  in  which  no 
sliding  was  permitted  between  the  specimen  and  the  steel  end  caps. 

The  simulated  specimen  is  1 1.43  cm  (4.5  in.)  in  length  with  a  diameter  of 
5.04  cm  (2  in.).  The  permeability  of  the  limestone  was  1 .03  x  lo’ m^.  which 
is  a  value  that  insures  a  uniform  pore  pressure  Held  throughout  the  mesh 
during  both  the  drairMd  and  undrained  simulations. 


SimuUitiort  of  drainod  triaxial  compraaaion  taat 

A  drained  triaxial  compression  test  at  a  confuting  stress  of  200  MPa  was 
simulated  in  the  following  manner.  First,  equal  incremenu  of  vertical  and 
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radiaJ  normal  stresses  were  applied  to  the  boundaries  of  the  grid  until  the 
stresses  equaled  200  MPa.  Second,  increments  of  vertical  stress  were  applied 
until  the  total  vertical  stress  reached  550  MPa.  Finally,  incremenu  of  venical 
stress  were  removed  until  a  hydrostatic  state  of  total  stress  was  obtained. 


The  output  from  this  calculation  is  ploned  in  the  form  of  contour  plots  of 
several  stress  or  strain  parameten,  i.e.,  >/}2d>  piutic  volume  strain,  axial 
strain  and  radial  strain  (Figures  6. 13  and  6. 14).  The  contour  plots  preseru  the 
state  of  stress  or  strain  in  the  specimen  at  the  time  of  peak  total  vertical  stress 
(Note:  the  end  cap  is  not  included  in  these  contour  plots).  With  the  exception 
of  the  upper  15  to  20 


Tabla  6.1. 

Laboratory  Calculatad  Straaa  and  Strain  Valtiaa 


percent  of  the  speci¬ 
men.  i.e.,  near  the 
interface  of  the 
specimen  and  end  cap, 
the  state  of  stress 
within  the  specimen  is 
relatively  uniform 
(Figure  6.13).  This  is 
also  true  of  the  plastic 

volume  strains  within  the  specimen  (Figure  6. 13).  However,  both  the  axial 
and  radial  strains  (Figure  6. 14)  exhibit  significant  gradiems  throughout  the 
specimen.  The  smallest  axial  strains  (0.03  m/m)  are  at  the  top  of  the 
specimen;  the  largest  (0.20  m/m)  develop  at  the  center  of  the  specimen.  The 


Axial  Stram 

IS.8% 

Radial  Svain 

•6.2% 

Princtpal  Straa*  Oiffaranca 

310  MPa 

>/ J20 

1 79  MPa 

•  • 
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radial  strains  vary  from  approximately  zero  at  the  specimen-end  cap  interface 
to  as  much  as  -0.065  m/m  at  the  center  of  the  specimen. 

Table  6. 1  contains  values  of  stress  and  strain  that  would  be  calculated  from 
laboratory  measured  load  aiKl  deformation  measurements.  Stresses,  e.g., 
principal  stress  difference,  were  corrected  for  the  changing  cross-sectional 
area  of  the  test  specimen.  The  stress  values  underestimate  the  strength  of  the 
test  specimen  by  ^proximately  5  percent,  179  MPa  (from  above  table)  versus 
190  MPa  (average  stress  throughout  specimen).  The  axial  strain  of 
15.8  percent  represents  only  a  small  portion  of  the  calculated  axial  strain 
within  the  test  specimen,  while  the  radial  strain  of  -6.2  percent  is 
representative  of  the  radial  strains  in  the  central  portion  of  the  test  specimen. 

This  calculation  implies  that  the  state  of  stress  within  the  test  specimen  is 
net  significantly  effected  by  end  c^  restraints.  Uniform  stresses  occur 
throughout  major  portions  of  the  specimen,  in  contrast,  large  axial  and  radial 
strain  gradients  were  developed  in  the  test  specimen.  This  implies  that  some 
type  of  end-cap  lubrication  should  be  used  if  uniform  states  of  strain  are 
desired. 


Simulation  of  conaolidatad  undrainad  triaxial  comprassion  taat 

A  consolidated  undrained  triaxial  compression  test  at  a  confining  stress  of 
200  MPa  was  simulated  with  the  FE  code  JAM.  To  begin  the  calculation, 
equal  increments  of  vertical  and  radial  normal  stresses  were  applied  to  the 
boundaries  of  the  grid  until  a  hydrostatic  stress  of  200  MPa  was  achieved. 
During  this  hydrostatic  loading,  pore  fluid  was  allowed  to  drain  from  the 
specimen.  The  boundary  conditions  were  then  changed  so  that  no  pore  fluid 
could  drain  from  the  specimen.  Finally,  increments  of  vertical  stress  were 
applied  until  the  solution  would  not  converge,  which  implied  that  the  specimen 
had  failed.  Failure  occurred  at  a  total  axial  strain  of  aj^roximately 
4.7  percent.  A  uniform  pore  fluid  pressure  existed  throughout  the  specimen. 

The  output  from  this  calculation  is  plotted  in  the  form  of  contour  plots  of 
\p2D'  volume  strain,  axial  strain,  and  radial  strain  (Figures  6.15  and  6.16). 
The  contour  plots  present  the  state  of  stress  or  strain  in  the  specimen  at  the 
time  of  specimen  failure  (Note:  The  end  cap  is  again  not  included  in  these 
contour  plots).  The  Table  6.2. 

^2D  contours  (Figure  Laboratory  Calculated  Strata  and  Strain  Values 
6.15)  illustrate  that  a  for  Undrainad  TXC  Test 
uniform  state  of  stress 
exists  within  a  majority 
of  the  test  specunen; 
significant  gradients 
only  exist  near  the 
specimen-end  cap  inter¬ 
face.  The  same  is  true 
of  the  volume  strain 


Akial  Strain 

4  7%  1 

Radial  Strain 

0.4%  1 

1  Vokima  Strain 
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contours.  The  calculated  volume  strains  indicate  that  the  specunen  was 
conqjacting.  The  calculated  pore  pressures  confirm  this  observation,  i.e.,  they 
increase  continuously  until  spec  ten  failure  occurs.  Due  to  the  small  axial 
strain  level  at  which  specimen  failure  occurs,  significant  gradients  of  axial  and 
radial  strain  were  not  developed  in  the  test  specimen  (Figure  6.16).  Axial 
strains  vary  h’om  0.02  m/m  at  the  top  of  the  specimen  to  less  than  O.OSS  m/m 
in  the  center  of  the  specimen.  Radial  strains  range  from,  approximately  zero 
at  the  specimen-end  cap  interface  to  -0.0C4  m/m  at  the  center  of  the  specimen. 

Table  6.2  contains  values  of  stress  and  strain  that  would  be  calculated  from 
laboratory  measured  load  and  deformation  measurements.  As  in  the  drained 
simulation,  stresses  were  corrected  for  the  changing  cross-sectional  area  of  the 
test  specimen.  The  laboratory  calculated  stress  values  correspond  with  the 
values  h’om  the  test  specimen  simulation,  i.e.,  124  MPa  (from  above  uble) 
versus  123  MPa  (average  stress  throughout  specimen).  The  laboratory  calcu¬ 
lated  volume  strain  of  3.8  percent  underestimates  the  simulated  volume  strains 
that  vary  between  4  and  4.4  percem  throughout  most  of  the  test  specimen. 

The  axial  strain  of  4.7  percent  agrees  with  the  calculated  axial  strain 
throughout  a  major  portion  of  the  test  specun«i.  The  radial  strain  of  -0.4%  is 
representative  of  the  simulated  radial  strains  in  the  central  portion  of  the  test 
specimen. 

This  calculation  demonstrates  that  significant  stress  and  strain  gradients  are 
not  developed  in  the  limestone  when  the  specimen  fails  at  small  axial  strains, 
despite  the  introduction  of  end  cap  restraint,  in  addition,  the  stresses  and 
strains  calculated  from  laboratory  measurements  correlate  well  with  the  actual 
stress  and  strain  sutes  within  the  test  specimen. 


Summary 

Numerical  simulations  of  limestone  behavior  under  drained  and  undrained 
boundary  conditions  were  presemed  in  this  chapter.  The  ability  of  the 
14-paraineter  cap  model  to  simulate  the  drained  behavior  of  the  limestone  was 
demonstrated  by  comparing  calculued  responses  of  hydrostatic,  uniaxial 
strain,  and  triaxial  compression  loadings  with  measured  or  recommended 
limestone  responses.  The  ability  of  the  FE  code  to  calculate  the  undrained 
behavior  of  the  limestone  was  demonstrated  by  comparing  calculated  respon¬ 
ses  of  uniaxial  strain  loadings  with  recommended  limestone  responses. 

Finally,  both  drained  and  undrained  TXC  test  siimdaiions  were  documented 
which  demonstrate  the  utility  of  the  FE  code  in  analyzing  laboratory  test 
sp(x;imen  conditions. 
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Summary 


This  rqwrt  (^<ocunients  the  features  and  algorithms  implemented  into  the  FE 
code  JAM.  The  FE  code  JAM  is  a  numerical  tool  with  the  capability  to: 

•  calculate  strains,  total  and  effective  stresses,  and  pore  fluid  pressures 
for  fully-  and  partially-saturated  porous  media. 

•  calculate  the  time  dependent  flow  of  pore  fluids  in  porous  media. 

•  model  nonlinear  irreversible  stress-strain  bdiavior.  including  cotqiled 
shear-induced  volume  change,  and 

•  simulate  the  effect  of  nonlinear  pore  fluid  compressibility  and  the 
contribution  of  the  compreuibility  of  the  grain  solids  for  stresses  up  to 
several  hundred  megapascals. 

in  this  r^rt.  the  FE  model  in^lemented  into  JAM  was  described,  and 
equations  were  developed  for  the  residual  forces.  The  features  of  the  cap 
model  and  the  relevent  etpiatkms  were  documented,  and  the  steps  required  to 
implement  the  cap  model  into  the  FE  code  were  summarized.  Other  consti¬ 
tutive  models  available  in  the  code  were  also  reviewed. 

The  equations  of  Mate  for  air,  water,  and  the  grain  solids  were  documented, 
and  the  equations  for  the  compressibility  of  an  air-wuer  mixture  were  devel¬ 
oped.  Several  documented  veriflcation  problems  demonstrated  that  the  pro- 
grim  works  correctly.  These  problems  included  one-  and  two-dimensional 
ctmsolidatioo  problems.  Cryer's  problem  of  a  consolidating  q)here  of  soil,  and 
a  thick-walled  cylinder  problem. 

Numerical  simulations  of  limestone  befaavim  under  drained  and  undrained 
boimdary  conditioni  were  presented.  A  14-parameter  cap  model  rooddled  die 
skeletal  properties  of  the  limestone.  Single  elemeni  calculations  demonstrated 
the  ability  of  the  FE  code  to  simulatt  both  the  drained  and  undrained  re¬ 
sponses  of  the  limettone  under  several  different  load  and  unload  boundary 
conditions.  The  utilhy  of  the  FE  code  was  donomtrated  by  the  simulaikm  of 
drained  and  undrained  triaxial  compressioo  teso.  and  the  i^uence  end  c^> 
rmtraini  had  on  the  stress  and  nrain  stales  in  the  test  specimen. 
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